- Email: [email protected]
Discrete linearized least-squares rational approximation on the unit circle
__ __ l!.F!
CQa ELSEVIER
JOURNAL OF COMPUTATIONAL AND APPLIED MATHEMATICS
Journal of Computational and Applied Mathematics 50 ( 1994) 545-563
Discrete linearized least-squares rational approximation on the unit circle Marc Van Bare1 *, Adhemar Bultheel Department
of Computer Science, Katholieke
Universiteit Leuven, Celestijnenlaan
2OOA, B-3001 Heverlee, Belgium
Received 10 August 1992; revised 5 February 1993
Abstract We already generalized the Rutishauser-Gragg-Harrod-Reichel algorithm for discrete least-squares polynomial approximation on the real axis to the rational case. In this paper, a new method for discrete least-squares linearized rational approximation on the unit circle is presented. It generalizes the algorithm of Reichel-AmmarGragg for discrete least-squares polynomial approximation on the unit circle to the rational case. The algorithm is fast in the sense that it requires order ma computation time where m is the number of data points and LY is the degree of the approximant. We describe how this algorithm can be implemented in parallel. Examples illustrate the numerical behavior of the algorithm. Keywords: Rational interpolation; Rational approximation; Linearized least squares; (Block) unitary Hessenberg; (Block) Schur parameters; (Block) orthonormal polynomials; Parallel
1. Introduction In our previous publications [ 14,151, we investigated the rational interpolation problem not only for scalar function values but also for the vector case. The standpoint we took was a theoretical one. We found a parametrization of all solutions of the rational interpolation problem and fast algorithms to compute this parametrization. Fast means of order m2 floating-point operations (flops) if there are m interpolation data. However, the algorithms which we designed are not numerically stable. This limits their applicability to theoretical investigation of the rational interpolation problem or to those cases where there are no round-off errors, e.g. when the computations are done within a finite field. In our paper [ 161, we replaced the interpolation conditions by a discrete least-squares criterion and described an algorithm to solve this problem when all interpolation points are taken on the real axis. It generalizes the algorithm of Reichel [9] for the polynomial case. The latter is inspired * Corresponding author. e-mail: [email protected]. 0377-0427/94/$07.00
@ 1994 Elsevier Science B.V. All rights reserved
SsDr0377-0427(93)E0160-N
546
M. Van Barel, A. BultheeWJournal
of Computational and Applied Mathematics
50 (1994) 545-563
by the Rutishauser-Gragg-Harrod algorithm [ l&12] for the computation of Jacobi matrices. A similar way of reasoning was followed in [ 4,6]. In this paper, we construct an algorithm solving the complementary problem, i.e., when all interpolation points are taken on the unit circle. This algorithm generalizes the algorithm given by Reichel, Ammar and Gragg [lo]. The latter is inspired by the inverse unitary QR algorithm for computing unitary Hessenberg matrices [ 21. It will turn out that our algorithm is fast, more specifically it is of order m(~ flops where (Y indicates the degree of the rational approximant. Implementing this algorithm in a pipelined fashion on (Y processors working in parallel reduces the computational time to order m.
2. The theory behind the algorithm Let us assume that the function values are given in the form filei for the abscissae zi, i = 1,2,..., m, where fi, ei and Zi are all complex numbers. If we define the degree of a rational form n(z)/d(z) with (n,d) Z (0,O) as deg(n(z),d(z)) = max{degn(z),degd(z)} (deg0 = -l), then in the rational interpolation problem one wants to describe all polynomial couples (n( z ) , d( z ) ) of minimal degree which satisfy the interpolation conditions n(zi)
fi
i =
1
d(~i)=z,
2
,
,...,m.
In [ 141, a parametrization was given and an efficient algorithm to solve this problem. When the data are corrupted with noise, one does not want the interpolation conditions to be satisfied exactly. Definition 2.1 (Proper rational approximation problem). Given the data points zi with corresponding estimated function values filei, i = 1,2, . . . , m, given the degree cy and the weights ,jp) > 0, i = 1,2,..., m, we look for a rational form n( z )/d( z ) of degree < cx satisfying the following least-squares approximation criterion: minimize
distFpj (n, d) = 2
$p’ 1rip’ j2,
i=l
with ,!” = fi/ei - n( zi)/d(zi), ]r!‘)]* = rjp)ri(p) and where dist(,, (n, d) denotes between the rational function n( z ) /d( z ) and the data.
the l2 distance
The proper rational approximation problem is a nonlinear least-squares problem which can only be solved in an iterative way. Therefore, we rather minimize the norm of the linearized residual vector with components ri =
fid(Zi)
-
ein(Zi),
i= 1,2 ,...,
m.
We shall fix the degree of the approximant n(z)/d( z) to be LY= max{deg n( z), deg d( z)} where usually (Y << m. We also normalize the approximant in the following sense. Suppose n(z) = no + and d(z) =do+diz +...+d,z”; then we require n, = 1 if degn(z) = (Y and nlz + ..- +n,z” d, = 1 if deg n( z ) < cx (and deg d (z ) = a). Note that we have given preference to the numerator polynomial in the normalization procedure. However, by switching the role of fi and ei and of n (z )
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
50 (1994) 545-563
and d( z ), we get the symmetric problem giving preference to d( z ) in the normalization Thus in this paper, we shall solve the following rational approximation problem.
547
procedure.
Definition 2.2 (Linearized rational approximation problem). Given the data points zi with corresponding estimated function values j’i/ei, i = 1,2, . . . , m, given the degree a and the weights wi > 0, m, we look for the normalized rational form n( z ) /d( z ) of degree cy satisfying the i = 1,2,..., following least-squares approximation criterion: minimize
dist2(,,d)
= kwilri]‘, i=l
with Ti = f,d(zi) form (n(z),d(z))
- ein(zi), lri/* = cri and where dist(n, d) denotes the distance between the rational and the data.
Note that the function values fi/ei can be replaced by kifi/( kiei) (with ki # 0). This yields a different value of the residual ri. Solving the linearized rational approximation problem with k,f,, k,ei instead of fi and e, is equivalent to solving the problem with the original fi, ei but with the weights Wi]lii]* instead of Wi. Choosing ki, i = 1,2, . . . , m, can be done in several ways, e.g., one ,..., m. could normalize the data such that lfij2 + leil* = 1 orsuchthatmax{lf,),lei]}=1,i=1,2 The solution of the linearized problem can be used to obtain a solution of the proper problem as follows. Suppose we know the values d(P) (zi)) i = 1,2, . . . , m, with rz(P)/d(P’ a solution of the proper problem. If we solve the linearized problem with weights
we get (n(P), d(P)). However, in practice, we do not know the values d’P’( zi). In this case we can make an estimation of these values, compute the solution of the linearized problem, take the denominator of this solution as a new estimation of the final d(p), and so on. This algorithm was proposed by Loeb for the 1, norm and by Wittmeyer for the l2 norm [ 31. In Examples 7.2 and 7.3, we shall show the influence of executing one iteration step of this algorithm. Of course one could also use the solution of the linearized problem as a starting value for other iterative schemes. Let us now look how the linearized problem can be solved in an efficient way. The linear leastsquares problem can be reformulated as a block problem, i.e., we should minimize
with
Denoting
the approximant
as a( z ) = [ d( z > n( z ) IT, this can be formulated as
M. Van Barel, A. Bultheel/Journal
548
of Computational and Applied Mathematics 50 (1994) 545-563
with
Il4z)ll: = (4Z)dZ))”
and
(a(z),b(z)),=
emTuib(Zi), i=l
for any a( z ) , and b( z ) polynomial A(z),B(z) E (C[Z]~‘~, we defi?re (A(z),B(z))u
vectors
in C [ z ] 2. We also need a matrix
version.
So, if
= eA(zi)‘uiB(zi)* i=l
Now suppose that we can construct some orthonormal polynomials { &,, c#J~,. . .} in C [ z ] 2x2, such . + &, with &k regular and upper triangular and such that that&(z) =&zk+e. (~i(z)~~j(z))u
=6ijz2v
with Z2 the 2 x 2 unit matrix. Then we can express a( z ) as a linear combination
a(z)
Ui E C2.
=e+i(Z)Uip i=O
Thus,
k&O
\
/
i=l
Because a, cannot be zero,butallak,k=O,l the minimizing a( z ) ,
k,l=O
,...,
k=O
k=O
LY- 1, can be chosen as zero, we find that for
a &z)
=
z+k(Z)ak
=
+a(z)a,,
where a, is chosen to minimize ]]a(~) 11:= aTa,
= ]pu12 + ]r,12,
with a, = [pa rUIT.
Given the normalizing conditions for the approximant and the leading coefficient of tin (z ), this is an easy problem to solve. Our first task will be to construct the orthogonal block polynomials @,i. We use a generalization of the method introduced in [9,10] to compute a polynomial approximation in the least-squares sense.
M. Van Barel, A. B&heel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
549
3. The algorithm The algorithm starts by applying some unitary similarity transformations arranged in an array as follows:
on the data which are first
where we used the notation
fI
fi! . . * f:,
-e: -ei . . . _eo 1
T ’
Z=diag(zl,z2,...,z,,).
We have a double column in the SW quadrant, which is reflected as complex conjugate transposed rows in the NE quadrant, and in the SE quadrant we originally find the diagonal matrix of abscissae. We now apply unitary similarity transformations to the left and right of this initial matrix M to transform it into Q’HMQ’, a block upper Hessenberg matrix (the blocks are 2 x 2)) where Q’ = Z2$ Q. Therefore, if we take all abscissae zi on the unit circle, i.e., Izi( = 1, i = 1,2, . . . , m, QHZQ will be unitary (block upper Hessenberg) . When the data points lie on the unit circle, Reichel, Ammar and Gragg [ lo] solve the least-squares polynomial approximation problem using a Schur parametrization of the unitary Hessenberg matrix that plays the central role in the theory. Let us return to the computational scheme and explain how the transformation Q’HMQ’ is actually performed using order m* computational work. Therefore, we need the following properties of unitary block Hessenberg matrices. Theorem 3.1 (Block Schur parametrization). m can be represented as H=G,G2G3...
Every unitary block Hessenberg matrix H of dimension
Go,,- I 3
with m’ = 1V$(m + 1) J, the greatest natural number smaller than or equal to i (m + 1). Each Gk is a unitary matrix of the form
r 12(k-1)
Note that Ik is the k x k unit matrix. For k < 0, this represents the empty matrix. The matrices yk, (7k, jjk and &k (block Schur parameters) are 2 x 2 except when m is odd and k = m’ - 1. In this case y~,~_lis2x2,a,~_lis1x2,6,~_~is2x1and~,~_~is1x1.Thematricesak,k=1,2,...,m’-1, are the subdiagonal blocks of the unitary block Hessenberg matrix H. Proof. Each unitary block Hessenberg H= [
-Y1
x
ul
x
0
x
1 .
matrix H can be written as
M. Van Barel, A. Bultheel/Journal
550
of Computational
and Applied Mathematics
50 (1994) 545-563
Because HHH = Z,,, we have that yyy~ + aycr, = Z,. Hence, there exist matrices 9, and &, such that
G,=
[2
;:
is unitary. Therefore, GyH=
[
;
;,
J
GYH is unitary and equal to
I
.
Because HHH = I,,,, XXH = 0 or X = 0. Hence, we can continue the same way of reasoning considering the smaller unitary block Hessenberg matrix H’ to get G2 and so on. Note that the subdiagonal blocks of H’ are the same as those of H skipping the first one. 0 In the sequel, we also need the following theorem. Theorem 3.2. Given the 1 x 2 matrix A and the 2 x 2 matrix B, there always exists a unitary matrix T such that
[a] =[OZJ
l-H
9
with C a 2 x 2 upper triangular matrix. Proof. A possible choice for T follows by taking the QR factorization
of the matrix [A” BHIH.
0
We could transform the initial array M into a unitarily similar block Hessenberg matrix using Givens reflections (or rotations) or Householder transformations. However, this would require 0( m3) computational work. Because H = QHZQ is unitary block Hessenberg, we can parametrize this matrix using block Schur parameters. We now construct a recursive algorithm where we transform these block Schur parameters in each step adding one new data point A = zmil, fm+i, emfl, w,+~. After m steps, the initial data matrix M is transformed into
OCT,” 0
Q'"MQ'=
H-1 (TO
H
,
0
with H = QHZQ a unitary block Hessenberg matrix with upper triangular subdiagonal (as well as ao). In step m + 1, we add the (m + 1)st data point . ., (T,,_~
~I,~Z,.
with ?ro = [J-L+, - ek,,, ] which is unitarily similar to
blocks
M. Van Barel, A. BultheeUJournal
Using Theorem
50 (1994) 545-563
551
such that Tt
with c(, upper triangular.
,p_J
[:
Using a block Schur parametrization
L]
= [7
ZJL]
which is again unitary block Hessenberg.
[7
and Applied Mathematics
3.2, we derive a unitary matrix TO:
To =
[Y
of Computational
ZnF*] [:
:I
[?
[6’ r;l%‘[A Hence,
for H = GIG2 . . . G,,,_,, we can write
FJ
[” ,,I is unitarily similar to
ZL]?
which is unitary block Hessenberg
except for the principal 3 x 3 block submatrix.
With the notation
[j~$j=[-$j$][6\J, we can compute this 3 x 3 block or 5 x 5 principal submatrix as
[”
$
,]
[’
;;I
;;I
[“‘i
Zz]
#PO - P$l vo -a;fio - PoHYlPO p,“*r &$, -&.L() - B,Hy,vo po”# - #y*&
= [
UIVO
71
1
-Pul fi, VI81 1
Using Theorem T, =
UlfQO
.
3.2, we can always find a unitary matrix T,: such that Tp
[
with a{ upper triangular and 7rl = -&‘,uo - &$qvo. Multiplying the 5 x 5 principal submatrix to the left by Z2 @ T,“, we get a unitary block Hessenberg matrix which can be written with -7; = CY,HZ_L~ - #y, v. and the g{ from above using block Schur parameters as -y;
5;
4 ?:
I To get a unitarily
12
1 Ii
-$ PF1 P,
Gr 1
similar matrix, we have to multiply the unitary block Hessenberg
matrix
M. Van Barel, A. Bultheel/Journal of Computational and Applied Mathematics 50 (1994) 545-563
552
to the right by I2 CDT, 63 Zm_4. Hence, we can repeat the same reasoning for Tl, T2, T3, . . . as for To. Before the final step of the algorithm is executed, we have the following unitarily similar (m + 1) x (m + 1) matrix: 12(m’-2)
12(r?*‘-2) -afit,-
G;G; . . .G’,,_,
-rum
P,“L,
8;tL2
-2
Vml-2
t&2
cd
-2
j&-2
’ LS I
A5
with S = 1 if m is odd, and 8 = 2 if m is even. The unitary matrices GL are built up using the prime block Schur parameters. The final step transforms the last (3 + 8) x (3 + 6) block
-+$-2 &‘-2
B =
$1:
,1
[’
;?:’
;;J
[7:;’
z:
z8].
[ l
If m is odd, there exists a unitary 2 x 2 matrix T,,r_l such that
with q,,,_, = -&_2,um~_2
- i;trf,_2ymj--l v,,,/-2 and al,,__, upper triangular. Finally we get
The last unitary matrix Gk,_, is built up by these block Schur parameters. l If m is even, there exists a 3 x 3 unitary matrix T,r_, such that TEL
[ g2;;,_2]
= [+I
)
with r,,,> _, as above and c$,, _, upper triangular. Finally we get
[
z2T;,_,] B [”
Tm,_l] = [ iFA’
These block Schur parameters
4. The solution
$1:
l]
[ z2 ;t
;v]
1
give us the final unitary matrices GL,_, and GL?,.
of the rational approximation
problem
After applying m steps of the algorithm constructed [u‘ OIT with H unitary block Hessenberg
above, we get Q”ZQ = H and Q” [f’
given by its block Schur parameters
- e’] =
yk, (Tk, &_kand Tk,
M. Van Barel, A. Bultheel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
k= 1,2,. . . , m’ - 1. For the moment we also assume that all gk are nonsingular, the regular case. If some flk is singular, we call it the singular case.
553
which we shall call
4.1. The regular case When the block Schur parameters gk are regular, the block Schur parameters allow us to define the orthonormal block polynomials +k in a recursive way. First, we need the following properties for the block columns of Q based on the same Schur parameters. Theorem 4.1. The block columns of Q satisfy the recursion ti!k+lck
=
ZQk
with initialization
+
Q;+&
&7’k,
Ql = Qi = [f’
= zQk7’: + Q;,
k=
1,2,3 ,...,
m’-
1,
- e’] a;‘.
Proof. From Q” [f’ - e’] = [(T: OIT, we get that Q, = [f’ - e’] ~0’. We set Q{ = Qt. From Q”ZQ = H, we get ZQ = QH = QGIG~.. . G,,! _, . Denote the (k + 1) st block column QGl G2Gk by Q:,, . Hence, by comparing the block columns of ZQ = QH, we get the recursion
of
Q;,, = Q:h + Qk+l%.
ZQk = -Q:?‘k + Qk+tak,
We get Qk+l out of the first equation and use this result to eliminate Qk+r from the second equation. Hence,
Q:,, = ZQkg;‘?k + Q;@-, +
Qk+l = (ZQk + Q;Yk)6’,
Ykg;‘?kh
Because [;T
;j
[;$
we can rewrite ~k’fk
$]=
[‘a I”i],
as $&L”
and &k + rkg;‘&
as 3,“.
This proves the theorem.
We define the block polynomials &, k = 0, 1,2, . . . , by a similar recursion block polynomials are orthonormal with respect to the scalar product (A(z),
B(z)),, =
•!
and show that these
~~TsB(z,,. i=l
Definition 4.2. We define the block polynomials based on the block Schur parameters: 4kcZ)u.k
= i$k-I
(z>
+
&&)Yk,
~~(z)~~=Z~k-I(Z)Y~+~:-,(Z)r
with the initialization Property 4.3.
C$kand c&, k = 0, 1,2, . . . , by the following recursion
k= 1,2 ,...,
m’-
(4.1)
1,
& = &, = u;‘.
The block polynomials
+k and 4; satisfy the following
properties
(with &,, = a:).
554
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
50 (1994) 545-563
l The polynomial qf~k has degree k, highest-degree coejficient u;‘u;’ ’ ’ ’ ~;?,a; and constant coefficient &;“LY;” +. . &;T, &F”@. a The polynomial q5: has degree k, highest-degree coejficient a;‘~;’ e. . a;‘Tk and constant coefj’icient a;“&,” . . . &iT, if?;“. l The block polynomials & and 4; are connected to the block columns Qk and Qi as follows:
Q;,, = D@s:,
Qk+l = D@kr with D = diag{ [fi jined similarly).
- e{]
, . . . , [f:,
-e;,]}
and
@k
=
[+k(zdT
4k(zdT
. . . +k(zrdTIT(@;
is
Proof. The properties immediately follow from the recurrence relations for the block polynomials and 4: and the fact that Yk, (7k, &_k,Tk are block Schur parameters. 0 The following result allows us to solve the discrete linearized problem on the unit circle. Theorem 4.4 (Orthonormality). The block polynomials respect to the scalar product (., *)u, i.e., (4i(z), @j(z)),
least-squares
&
rational approximation
C$kas dejined above are orthonormal =
de-
with
SijZz.
Proof. Using the results of the previous theorem, the scalar product can also be written as (41 (Z) 9rb;(Z >)c = @iDHD@; = Q,“,,Qj+I = SijZIv because QHQ = I,,,.
0
Let us now return to the rational approximation problem. To solve this problem, we shall write the solution in terms of the orthonormal block polynomials 4j. Because the highest-degree coefficients of these block polynomials are upper triangular, we can always write our solution n( z ) /d (z ) of degree LYas follows:
with (p,, 7,) # (0,O). Remember that the solution we are looking for has to be normalized. Property 4.3, we know that the highest-degree coefficient of 4a is equal to
From
(4.2) The specific value of x is not important for the moment. To normalize n( z ) /d( z ), we have to choose 7, = n,“, S2.j if n, # 0 or per = nj”, Sl,j if n, = 0. To solve the rational approximation problem, we have to minimize
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
SO (1994) 545-563
555
As we have seen in Section 2, we can make this entity smaller by putting &. and r/, equal to zero for (n(z),d(z)) andthedatato (lp,]*+]r,]*)‘/*. k=O, 1,. ..,a-l.Thisreducesthedistancebetween The way in which we have stated the rational approximation problem requires n(z ) /d( z ) to have degree cy. From the normalization conditions, it is clear that we should therefore take the “right” structure, i.e., 7, = I-I S2.i and pa = 0,
if IJ
and the “left” structure,
Pa =fi
Is2,jl6 I-J ISl.jI9 j=O
j=O
j=O
i.e.,
Si,j and 7, = 0,
if fi
IS*,jl3 fi
ISl,jl,
j=O
.i=O
.i=o
where, in case of equality of the two distances, we can choose between the two solutions. Note that this choice corresponds to our problem setting, but that we have the possibility to solve also slightly different problems in an equally simple way. For instance, we could have required n( z ) /d( z ) to be strictly proper, i.e., with IZ, = 0 and d, = 1. Then we had to choose the second solution, i.e., the one with 7, = 0. It is easy to generalize this to more complicated structures than the strictly proper one. Suppose for example that we look for a solution with ~1, = 0, IZ,_~ = 0 and d, = 1. From n, = 0, we easily derive that 7, = 0. From d, = 1, it follows that pa = nj”, sl,j. To make n,_i equal to zero, we have to take T,_~ equal to a specific value, possibly different from zero. We shall not follow this path in the sequel but concentrate our attention on the original problem setting. The other cases are quite similar. Let us investigate how good the rational approximation is. in terms of the parameters pk, ?-k, Theorem 4.5. If we write our rational approximant n(z)/d(z) k = 0,1,2 ,..., cy, then the approximation error made in the ith point z;, i = 1,2, . . . , m, can be measured as follows:
f,‘d(zi> - e(n(zi>= ~QG+I
[
:I] .
!A
Proof. This result follows immediately [j-i’ - e;]
@kc&)
=
QG+I.
the ith row of the block column Qk+i. For eXampk, if we take the Pn = n,“, Si,j. Th e approximation f:d(zi)
from the fact that
- eInCzi> = Qi,a+l
0
proper SOhtiOU, all pk and rk VaheS error in each point zi can thus be written as
Strktly
[
PU 0
I ’
i.e., the error is given by a scaled version of the first column of Qn+, .
are
zero
except
one,
M. Van Barel, A. Bultheel/Journal
556
of Computational and Applied Mathematics 50 (1994) 545-563
4.2. The singular case Until now, we have taken for granted that all the subdiagonal blocks uj were nonsingular. What will happen if one of the upper triangular matrices aj is singular, i.e., if si,j or s2,j is equal to zero?
Theorem 4.6. We shall assume that ga is the first block Schur parameter s2,a equal to zero. We de&e the block polynomial & as follows: #a = IfSl,a
=
Z&l
(2
> +
4&_,
(z
with ~1,~ or (if it exists)
ha.
0, then
gives us a normalized rational interpolant Otherwise, ~2,~ = 0 and
n( z ) /d( z ) with d, = 1.
gives us a normalized
n( z ) /d (z ) with n, = 1.
rational interpolant
Proof. Because the previous block Schur parameters compute
c$~ based on $a-1
and &_, . Moreover,
aj, j = 0, 1,2, . . . , a - 1, are regular, we can c$~ has degree (Y with the same highest-degree
coefficient as c$~_, , i.e., n,“;;’ fl,:‘. Using the relationship between the block columns of the unitary matrix Q and the orthonormal block polynomials @jy j = 0, 1, . . . , a - 1, we get
Qi,a+l u,=
[fi
-ejfja(Zi),
i=1,2
,...,
m,
with CT, =
[
Sl,o
S3,a
0
S2.a
If we summarize
1 *
the two cases as
we get that [fi
- ei]
tlz,’ [ I 1 =[fi’
- 414aC-G)
f
[I
=Qi,a+lua [f]
=O,
i=1,2
,...,
m,
for the specific choices of jj and ?. Hence, n( z)/d( z ) is a rational interpolant for the given data. By looking at the highest-degree coefficient, we see that this rational interpolant is normalized. This proves the theorem. 0
M. Van Barel, A. BultheeUJournal
of Computational
and Applied
Mathematics
50 (1994)
545-563
551
The rational function n( z ) /d( z ) from the previous theorem approximates the given data exactly, i.e., n( z)/d( z) is an interpolant of the given data. It is clear that this situation will appear rarely in practice when the data are only given with a certain finite precision.
5. Deriving a rational approximant
with real coefficients
In a lot of practical applications, e.g., digital filter design, model reduction, modal analysis, system identification, one looks for a rational approximant n( z ) /d( z ) with real coefficients. The algorithm, as described before, results in complex coefficients in general. If we want n( z)/d( z ) to have real coefficients, the data should be consistent with this restriction, i.e., if z is a given abscissa with corresponding estimated function value f/e, then f/a can be considered as an estimated function value in F. Using this fact, the initial data matrix M can be transformed using a simple unitary matrix into a similar real matrix. From this point on, all computations can be done using real arithmetic, resulting in real block Schur parameters, real block orthogonal polynomials and finally a real rational approximant. The simple unitary matrix to do this initial similarity transformation can be chosen for each data pair in several ways, e.g.,
because
6. Parallel implementation Regarding the amount of computational work, we can make a distinction between the computation of the unitary block Hessenberg matrix QHZQ from the data and the computation of the orthonormal block polynomials. Fixing the degree LY< m of the rational approximant implies that we can obtain the rational approximant from +a (z ) . Hence, it is not necessary to compute the complete block Schur parametrization but only the first part consisting of Gr , G2, . . . , G,. Compared to the complexity rn2 of computing the complete block Schur parametrization, this decreases the amount of computational work to ICY. The second part of the solution method requires the computation of &,+,, . . . , $J,, which requires a number of floating-point operations of the order of (Ye. Therefore, since (Y < m, it is important to decrease the execution time of the first part by parallel implementation. This could be done by a sort of pipeline method. Let us assume that we have cv < m processors. Each processor is coupled to a Gk, k = 1,2, . . . , a, matrix and updates the corresponding block Schur parameters in each step. Once the kth processor
M. Van Barel, A. Bultheel/Journal of Computational and Applied Mathematics 50 (1994) 54.5-563
558
IO-4
10-7 -
IO-10
IO-13IO-16
0
+ 10
5
1
IS
20
25
30
35
40
Fig. 1. The modulus of the function values of Example 7.1.
has computed CQ, Pk, &, & and Tk, the (k + 1)th processor can use this information to update the corresponding block Schur parameters. Once the pipeline is started up, the execution time is equivalent to O(m) floating-point operations.
7. Examples In this section, we give three examples. The first one is rather theoretical, the others are more practical. In the three examples, we look for a rational approximant with real coefficients (see Section 5). To measure the goodness of the rational approximant, we consider a relative and absolute approximation error in each point zi. We take all the weights wi(p) = 1. The absolute error that is plotted in the figures is then nothing else but Jrip)j. Similarly, when the relative error is plotted, we plot jr~P’/(fi/ej)l. The data is normalized in such a way that Cyl, (Ifil’ + lei12) = jlaoll; = 1 (11. IIF stands for the Frobenius norm). We can do this by computing Wi such that wi(lfi12
+
lei12)= $.
(7.1)
The value of the rational approximant we could also use the recurrence (4.1) the unit circle, as in the examples here, precision on a DEC workstation using
n( zi) /d( zi) is computed using the Homer scheme. Of course, where z is replaced by zi. For equally distributed points zi on one could also use FIT. All computations were done in double PRO-MATLAB with machine epsilon eps = 2.2204. 10-16.
Table 1
0
1 2 3 4 5
7.4349607 7.265 425 3 6.1803399. 6.1803399. 6.083 663 2 5.8550486.
’ 10-l . 10-l 10-l 10-l . 10-l 10-l
6.687403 1 5.5589297. 5.257311 1 6.305 926 4 5.473 828 7 . 5.123981 1.
10-l 10-l 10-l lo-l6 lo-l6 lo-l6
M. Van Barel, A. Bultheel/Journal
of Computational
and Applied Mathematics
50 (1994) 545-563
559
Table 2 Degree
Numerator coefficients
Denominator
0
3.152963 175 230421 IO-‘” 1.000000000000001 . loo 8.088600263652522. 10-l’ 1.000 000 000 000 000 10”
1.ooo 000 000 000 000 loo 1.633 234097 724558 . lo--l6 1.576481587615210. lo-l6 1.435 322 868 745 093 1O-l6
1 2 3
coefficients
Example 7.1. In the first example, we consider the m = 40 interpolation m, with function values i= 1,2,..., fi
= Z” + ’
q
Z.
points zi = exp(ij27r/m),
”
i.e., the polynomial z (Z - j) (z + j) evaluated at the points zi. The modulus of these function values is plotted in Fig. 1. The errors connected to the “left” and “right” solution of the linearized problem are given in Table 1. From this table, it is clear that we should choose cx equal to 3 and take the “right” structure for the rational approximant. Table 2 shows the coefficients of the numerator and denominator polynomial of the rational approximant of degree LY= 3. The absolute and relative errors are indicated in Fig. 2. We can say that the rational approximant we obtained here is very good, which, for this example, of course will not come as a surprise. In the next example circle.
we try to find a good approximant
for the exponential
function
on the unit
Example 7.2. The points zi are the same as in the first example. The function values are fi/ei = eXp(Zi), i = 1,2,. , ., m. When looking at the errors of the solutions of the linearized problem, we choose to take a “right” solution of degree 3. The absolute and relative errors are given in Fig. 3. By executing one iteration step of the Loeb-Wittmeyer algorithm adapting the weights w, as described in Section 2, we get the errors of Fig. 4. Note that, using the normalization of (7.1)) the relative error is almost a constant over the set of points zi. However, doing one iteration step, it is the absolute error having this property. A similar phenomenon can also be observed in the next example. Moreover, the winding number of the absolute error function exp( z ) - n( z ) /d( z ) on the unit circle
IO--
* 0
lo-)2 5
IO
15
20
25
30
35
40
* 0
5
10
15
Fig. 2. The absolute and relative errors of Example 7.1.
2”
25
30
35
1
40
M. Van Barel, A. BultheeUJournal of Computational and Applied Mathematics 50 (1994) 545-563
560
Fig. 3. The absolute and relative errors of Example 7.2.
is equal to deg IZ+ deg d + 1. In several examples, the error function for a best Chebyshev rational approximant n//d’ closely approximates a perfect circle about the origin having a winding number deg n’ + deg d’ + 1 [ 131. Hence, the rational approximant obtained here after doing one iteration step will be a close candidate to a best Chebyshev rational approximant on the unit circle. A more practical example is the following one. Example 7.3. As in the previous examples, we assume that the points zi are equally distributed on the unit circle, i.e., we take zi = exp(ij27r/m), i = 1,2, . . . , m, with m = 200. We choose in a certain random way five zeros with modulus between 0 and 0.99. Similarly, we choose three poles with modulus between 0.90 and 0.99 and two poles with modulus between 0 and 0.90. With these poles and zeros, a rational function is constructed with real coefficients. Hence, we add the complex conjugates of the poles and the zeros. The data fi/ei is equal to this rational function evaluated in zi, but we have added a uniformly distributed absolute error of 10P4. The modulus of the function values is shown in Fig. 5. Examining the errors connected to the linearized problem, we choose for a “right” solution of a reduced degree 8. The absolute and relative errors are given in Fig. 6. Doing one iteration as described in Section 2, we obtain the absolute and relative errors given in Fig. 7. The
+++*+*+*+**++*++++++++**+*+++**+*-
IO-J _f.ff.f
+
IO-J
*
* f
.
f
*+ +**
IO.4 0
5
10
I5
20
2s
30
35
40
Li7rT-J 0
Fig. 4. The absolute and relative errors of Example 7.2.
+
*
M. Van Barel, A. Bultheel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
561
102
++
: +
t
t 10’
* +++
5 * *+
100
/
* ; *
i
*
t
+ t +* + t* f* *f
t ** f : c
w++++
Fig. 5. The modulus of the function values of Example 7.3.
left side of Fig. 8 shows the poles (indicated by *) and the zeros (indicated by o) of the original rational function of degree 10 in the complex plane. The right side shows the poles and zeros of the rational approximant after doing one iteration step. Note that the dominant poles and zeros, i.e., the ones closely to the unit circle, are well approximated. These three examples show that the rational approximant approximates the given data quite well in the proper rational sense. Moreover, applying only one iteration step of the algorithm of LoebWittmeyer already gives an absolute error whose magnitude does not change very much anymore over the set of points zi. Because the amount of computational work is small (O(ma) and, in parallel, O(m) floating-point operations), the designer of a rational approximant has a lot of freedom in choosing the degree, the iteration method (how are the weights adapted), the number of iteration steps, . . . . The errors connected to the “left” and “right” solutions of the linearized problem give a good indication for the designer to choose between the several possible structures for the rational approximant (degree, “left” or “right”) . In the context of discrete-time linear systems, we can interpret the rational approximation problem solved in this paper as a system identification problem. The data fi/ei can be considered as
I
IO-5 ’
0
20
40
60
80
100
120
140
160
180
200
0
20
I 40
60
x0
Fig. 6. The absolute and relative errors of Example 7.3.
100
I20
140
I60
180
200
M. Van Barel, A. BultheeUJournal of Computational and Applied Mathematics 50 (1994) 54.5-563
562
Fig. 7. The absolute and relative errors of Example 7.3.
measurements of the frequency response of a discrete-time system. We are now looking for a linear model (a rational function) of limited complexity (small degree) such that the measured frequency response and the frequency response of the model almost coincide (taking into account the error on the measurements). For an introduction to system identification one could look at, e.g., [7, Chapter 81 and [ 11, Chapter 61. The rational approximation problem can also be seen as a model reduction problem. In the last example, we start with a model of degree 10 and we approximate this by a model of degree 8 [ 5, Section 191.
8. Conclusion We have developed an efficient method to solve the discrete linearized least-squares rational approximation problem. We have shown that this method can be implemented in a pipelined way on a parallel architecture. The numerical examples show that the computed rational approximant is also quite good as a solution of the proper rational approximation problem. It becomes a very good proper rational approximant after only one iteration step of the algorithm of Loeb-Wittmeyer. Because the algorithm developed here is so efficient, it could be used as the kernel of a graphical
-1
-0.5
0
0.5
I
Fig. 8. The poles and zeros of the original and the approximating
rational function of Example 7.3.
M. Van Barel, A. Bultheel/Journal
of Computational
and Applied Mathematics
50 (1994) 54.5-563
563
interactive program to design rational approximants for large data sets. As domains of application, we mention digital filter design, modal analysis, model reduction, system identification, . . . . In future work, this algorithm will be generalized to other structures for the rational approximant, i.e., the maximum degree of numerator and denominator polynomial can be chosen freely. Also more than two data columns f’ and e’ will be allowed, leading to block orthonormal polynomials with the dimension of the blocks bigger than 2.
References ‘[l] G.S. Ammar and W.B. Gragg, O(n*) reduction algorithms for the construction of a band matrix from spectral data, SIAM _I. Matrix Anal. Appl. 12 ( 199 1) 426-43 1. [ 21 G.S. Ammar, W.B. Gragg and L. Reichel, Constructing a unitary Hessenberg matrix from spectral data, in: G.H. Golub and P. Van Dooren, Eds., Proc. NATO Advanced Study Institute on Numerical Linear Algebra, Digital Signal Processing and Parallel Algorithms, Leuven, 1988 (Springer, Berlin, 199 1) 385-395. [ 31 1. Barrodale and J.C. Mason, Two simple algorithms for discrete rational approximation, Math. Comp. 24 ( 112) (1970) 877-891. [4] D.L. Boley and G.H. Golub, A survey of matrix inverse eigenvalue problems, Inverse Problems 3 (1987) 595-622. [5] A. B&heel and M. Van Barel, Padt techniques for model reduction in linear system theory: a survey, J. Comput. Appl. Math. 14 (3) (1986) 401-438. Updating and downdating of orthogonal polynomials with data fitting [6] S. Elhay, G.H. Golub and J. Kautsky, applications, SIAM J. Matrix Anal. Appl. 12 (1991) 327-353. [7] P. Faurre and M. Depeyrot, Elements of System Theory (North-Holland, Amsterdam, 1977). [8] W.B. Gragg and W.J. Harrod, The numerically stable reconstruction of Jacobi matrices from spectral data, Numer. Math. 44 ( 1984) 317-335. [ 91 L. Reichel, Fast QR decomposition of Vandermonde-like matrices and polynomial least squares approximation, SIAM .I. Matrix Anal. Appl. 12 ( 1991) 552-564. [lo] L. Reichel, G.S. Ammar and W.B. Gragg, Discrete least squares approximation by trigonometric polynomials, Math. Comp. 57 ( 1991) 273-289. [ 111 W.J. Rugh, Mathematical Description of Linear Systems (Marcel Dekker, New York, 1975). [ 121 H. Rutishauser, On Jacobi rotation patterns, in: N.C. Metropolis, Ed., ExperimentalArithmetic, High Speed Computing and Mathematics, Proc. Sympos. Appl. Math. 15 (Amer. Mathematical Sot., Providence, RI, 1963) 219-239. [ 131 L.N. Trefethen, Near-circularity of the error curve in complex Chebyshev approximation, _/. Approx. Theory 31 ( 1981) 344-367. [ 141 M. Van Bare1 and A. Bultheel, A new approach to the rational interpolation problem, J. Comput. Appl. Math. 32 (l&2) (1990) 281-289. [ 151 M. Van Bare1 and A. Bultheel, A new approach to the rational interpolation problem: the vector case, J. Comput. Appl. Math. 33 (3) ( 1990) 33 l-346. [ 161 M. Van Bare1 and A. Bultheel, A parallel algorithm for discrete least squares rational approximation, Numer. Math. 63 (1992) 99-121.
CQa ELSEVIER
JOURNAL OF COMPUTATIONAL AND APPLIED MATHEMATICS
Journal of Computational and Applied Mathematics 50 ( 1994) 545-563
Discrete linearized least-squares rational approximation on the unit circle Marc Van Bare1 *, Adhemar Bultheel Department
of Computer Science, Katholieke
Universiteit Leuven, Celestijnenlaan
2OOA, B-3001 Heverlee, Belgium
Received 10 August 1992; revised 5 February 1993
Abstract We already generalized the Rutishauser-Gragg-Harrod-Reichel algorithm for discrete least-squares polynomial approximation on the real axis to the rational case. In this paper, a new method for discrete least-squares linearized rational approximation on the unit circle is presented. It generalizes the algorithm of Reichel-AmmarGragg for discrete least-squares polynomial approximation on the unit circle to the rational case. The algorithm is fast in the sense that it requires order ma computation time where m is the number of data points and LY is the degree of the approximant. We describe how this algorithm can be implemented in parallel. Examples illustrate the numerical behavior of the algorithm. Keywords: Rational interpolation; Rational approximation; Linearized least squares; (Block) unitary Hessenberg; (Block) Schur parameters; (Block) orthonormal polynomials; Parallel
1. Introduction In our previous publications [ 14,151, we investigated the rational interpolation problem not only for scalar function values but also for the vector case. The standpoint we took was a theoretical one. We found a parametrization of all solutions of the rational interpolation problem and fast algorithms to compute this parametrization. Fast means of order m2 floating-point operations (flops) if there are m interpolation data. However, the algorithms which we designed are not numerically stable. This limits their applicability to theoretical investigation of the rational interpolation problem or to those cases where there are no round-off errors, e.g. when the computations are done within a finite field. In our paper [ 161, we replaced the interpolation conditions by a discrete least-squares criterion and described an algorithm to solve this problem when all interpolation points are taken on the real axis. It generalizes the algorithm of Reichel [9] for the polynomial case. The latter is inspired * Corresponding author. e-mail: [email protected]. 0377-0427/94/$07.00
@ 1994 Elsevier Science B.V. All rights reserved
SsDr0377-0427(93)E0160-N
546
M. Van Barel, A. BultheeWJournal
of Computational and Applied Mathematics
50 (1994) 545-563
by the Rutishauser-Gragg-Harrod algorithm [ l&12] for the computation of Jacobi matrices. A similar way of reasoning was followed in [ 4,6]. In this paper, we construct an algorithm solving the complementary problem, i.e., when all interpolation points are taken on the unit circle. This algorithm generalizes the algorithm given by Reichel, Ammar and Gragg [lo]. The latter is inspired by the inverse unitary QR algorithm for computing unitary Hessenberg matrices [ 21. It will turn out that our algorithm is fast, more specifically it is of order m(~ flops where (Y indicates the degree of the rational approximant. Implementing this algorithm in a pipelined fashion on (Y processors working in parallel reduces the computational time to order m.
2. The theory behind the algorithm Let us assume that the function values are given in the form filei for the abscissae zi, i = 1,2,..., m, where fi, ei and Zi are all complex numbers. If we define the degree of a rational form n(z)/d(z) with (n,d) Z (0,O) as deg(n(z),d(z)) = max{degn(z),degd(z)} (deg0 = -l), then in the rational interpolation problem one wants to describe all polynomial couples (n( z ) , d( z ) ) of minimal degree which satisfy the interpolation conditions n(zi)
fi
i =
1
d(~i)=z,
2
,
,...,m.
In [ 141, a parametrization was given and an efficient algorithm to solve this problem. When the data are corrupted with noise, one does not want the interpolation conditions to be satisfied exactly. Definition 2.1 (Proper rational approximation problem). Given the data points zi with corresponding estimated function values filei, i = 1,2, . . . , m, given the degree cy and the weights ,jp) > 0, i = 1,2,..., m, we look for a rational form n( z )/d( z ) of degree < cx satisfying the following least-squares approximation criterion: minimize
distFpj (n, d) = 2
$p’ 1rip’ j2,
i=l
with ,!” = fi/ei - n( zi)/d(zi), ]r!‘)]* = rjp)ri(p) and where dist(,, (n, d) denotes between the rational function n( z ) /d( z ) and the data.
the l2 distance
The proper rational approximation problem is a nonlinear least-squares problem which can only be solved in an iterative way. Therefore, we rather minimize the norm of the linearized residual vector with components ri =
fid(Zi)
-
ein(Zi),
i= 1,2 ,...,
m.
We shall fix the degree of the approximant n(z)/d( z) to be LY= max{deg n( z), deg d( z)} where usually (Y << m. We also normalize the approximant in the following sense. Suppose n(z) = no + and d(z) =do+diz +...+d,z”; then we require n, = 1 if degn(z) = (Y and nlz + ..- +n,z” d, = 1 if deg n( z ) < cx (and deg d (z ) = a). Note that we have given preference to the numerator polynomial in the normalization procedure. However, by switching the role of fi and ei and of n (z )
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
50 (1994) 545-563
and d( z ), we get the symmetric problem giving preference to d( z ) in the normalization Thus in this paper, we shall solve the following rational approximation problem.
547
procedure.
Definition 2.2 (Linearized rational approximation problem). Given the data points zi with corresponding estimated function values j’i/ei, i = 1,2, . . . , m, given the degree a and the weights wi > 0, m, we look for the normalized rational form n( z ) /d( z ) of degree cy satisfying the i = 1,2,..., following least-squares approximation criterion: minimize
dist2(,,d)
= kwilri]‘, i=l
with Ti = f,d(zi) form (n(z),d(z))
- ein(zi), lri/* = cri and where dist(n, d) denotes the distance between the rational and the data.
Note that the function values fi/ei can be replaced by kifi/( kiei) (with ki # 0). This yields a different value of the residual ri. Solving the linearized rational approximation problem with k,f,, k,ei instead of fi and e, is equivalent to solving the problem with the original fi, ei but with the weights Wi]lii]* instead of Wi. Choosing ki, i = 1,2, . . . , m, can be done in several ways, e.g., one ,..., m. could normalize the data such that lfij2 + leil* = 1 orsuchthatmax{lf,),lei]}=1,i=1,2 The solution of the linearized problem can be used to obtain a solution of the proper problem as follows. Suppose we know the values d(P) (zi)) i = 1,2, . . . , m, with rz(P)/d(P’ a solution of the proper problem. If we solve the linearized problem with weights
we get (n(P), d(P)). However, in practice, we do not know the values d’P’( zi). In this case we can make an estimation of these values, compute the solution of the linearized problem, take the denominator of this solution as a new estimation of the final d(p), and so on. This algorithm was proposed by Loeb for the 1, norm and by Wittmeyer for the l2 norm [ 31. In Examples 7.2 and 7.3, we shall show the influence of executing one iteration step of this algorithm. Of course one could also use the solution of the linearized problem as a starting value for other iterative schemes. Let us now look how the linearized problem can be solved in an efficient way. The linear leastsquares problem can be reformulated as a block problem, i.e., we should minimize
with
Denoting
the approximant
as a( z ) = [ d( z > n( z ) IT, this can be formulated as
M. Van Barel, A. Bultheel/Journal
548
of Computational and Applied Mathematics 50 (1994) 545-563
with
Il4z)ll: = (4Z)dZ))”
and
(a(z),b(z)),=
emTuib(Zi), i=l
for any a( z ) , and b( z ) polynomial A(z),B(z) E (C[Z]~‘~, we defi?re (A(z),B(z))u
vectors
in C [ z ] 2. We also need a matrix
version.
So, if
= eA(zi)‘uiB(zi)* i=l
Now suppose that we can construct some orthonormal polynomials { &,, c#J~,. . .} in C [ z ] 2x2, such . + &, with &k regular and upper triangular and such that that&(z) =&zk+e. (~i(z)~~j(z))u
=6ijz2v
with Z2 the 2 x 2 unit matrix. Then we can express a( z ) as a linear combination
a(z)
Ui E C2.
=e+i(Z)Uip i=O
Thus,
k&O
\
/
i=l
Because a, cannot be zero,butallak,k=O,l the minimizing a( z ) ,
k,l=O
,...,
k=O
k=O
LY- 1, can be chosen as zero, we find that for
a &z)
=
z+k(Z)ak
=
+a(z)a,,
where a, is chosen to minimize ]]a(~) 11:= aTa,
= ]pu12 + ]r,12,
with a, = [pa rUIT.
Given the normalizing conditions for the approximant and the leading coefficient of tin (z ), this is an easy problem to solve. Our first task will be to construct the orthogonal block polynomials @,i. We use a generalization of the method introduced in [9,10] to compute a polynomial approximation in the least-squares sense.
M. Van Barel, A. B&heel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
549
3. The algorithm The algorithm starts by applying some unitary similarity transformations arranged in an array as follows:
on the data which are first
where we used the notation
fI
fi! . . * f:,
-e: -ei . . . _eo 1
T ’
Z=diag(zl,z2,...,z,,).
We have a double column in the SW quadrant, which is reflected as complex conjugate transposed rows in the NE quadrant, and in the SE quadrant we originally find the diagonal matrix of abscissae. We now apply unitary similarity transformations to the left and right of this initial matrix M to transform it into Q’HMQ’, a block upper Hessenberg matrix (the blocks are 2 x 2)) where Q’ = Z2$ Q. Therefore, if we take all abscissae zi on the unit circle, i.e., Izi( = 1, i = 1,2, . . . , m, QHZQ will be unitary (block upper Hessenberg) . When the data points lie on the unit circle, Reichel, Ammar and Gragg [ lo] solve the least-squares polynomial approximation problem using a Schur parametrization of the unitary Hessenberg matrix that plays the central role in the theory. Let us return to the computational scheme and explain how the transformation Q’HMQ’ is actually performed using order m* computational work. Therefore, we need the following properties of unitary block Hessenberg matrices. Theorem 3.1 (Block Schur parametrization). m can be represented as H=G,G2G3...
Every unitary block Hessenberg matrix H of dimension
Go,,- I 3
with m’ = 1V$(m + 1) J, the greatest natural number smaller than or equal to i (m + 1). Each Gk is a unitary matrix of the form
r 12(k-1)
Note that Ik is the k x k unit matrix. For k < 0, this represents the empty matrix. The matrices yk, (7k, jjk and &k (block Schur parameters) are 2 x 2 except when m is odd and k = m’ - 1. In this case y~,~_lis2x2,a,~_lis1x2,6,~_~is2x1and~,~_~is1x1.Thematricesak,k=1,2,...,m’-1, are the subdiagonal blocks of the unitary block Hessenberg matrix H. Proof. Each unitary block Hessenberg H= [
-Y1
x
ul
x
0
x
1 .
matrix H can be written as
M. Van Barel, A. Bultheel/Journal
550
of Computational
and Applied Mathematics
50 (1994) 545-563
Because HHH = Z,,, we have that yyy~ + aycr, = Z,. Hence, there exist matrices 9, and &, such that
G,=
[2
;:
is unitary. Therefore, GyH=
[
;
;,
J
GYH is unitary and equal to
I
.
Because HHH = I,,,, XXH = 0 or X = 0. Hence, we can continue the same way of reasoning considering the smaller unitary block Hessenberg matrix H’ to get G2 and so on. Note that the subdiagonal blocks of H’ are the same as those of H skipping the first one. 0 In the sequel, we also need the following theorem. Theorem 3.2. Given the 1 x 2 matrix A and the 2 x 2 matrix B, there always exists a unitary matrix T such that
[a] =[OZJ
l-H
9
with C a 2 x 2 upper triangular matrix. Proof. A possible choice for T follows by taking the QR factorization
of the matrix [A” BHIH.
0
We could transform the initial array M into a unitarily similar block Hessenberg matrix using Givens reflections (or rotations) or Householder transformations. However, this would require 0( m3) computational work. Because H = QHZQ is unitary block Hessenberg, we can parametrize this matrix using block Schur parameters. We now construct a recursive algorithm where we transform these block Schur parameters in each step adding one new data point A = zmil, fm+i, emfl, w,+~. After m steps, the initial data matrix M is transformed into
OCT,” 0
Q'"MQ'=
H-1 (TO
H
,
0
with H = QHZQ a unitary block Hessenberg matrix with upper triangular subdiagonal (as well as ao). In step m + 1, we add the (m + 1)st data point . ., (T,,_~
~I,~Z,.
with ?ro = [J-L+, - ek,,, ] which is unitarily similar to
blocks
M. Van Barel, A. BultheeUJournal
Using Theorem
50 (1994) 545-563
551
such that Tt
with c(, upper triangular.
,p_J
[:
Using a block Schur parametrization
L]
= [7
ZJL]
which is again unitary block Hessenberg.
[7
and Applied Mathematics
3.2, we derive a unitary matrix TO:
To =
[Y
of Computational
ZnF*] [:
:I
[?
[6’ r;l%‘[A Hence,
for H = GIG2 . . . G,,,_,, we can write
FJ
[” ,,I is unitarily similar to
ZL]?
which is unitary block Hessenberg
except for the principal 3 x 3 block submatrix.
With the notation
[j~$j=[-$j$][6\J, we can compute this 3 x 3 block or 5 x 5 principal submatrix as
[”
$
,]
[’
;;I
;;I
[“‘i
Zz]
#PO - P$l vo -a;fio - PoHYlPO p,“*r &$, -&.L() - B,Hy,vo po”# - #y*&
= [
UIVO
71
1
-Pul fi, VI81 1
Using Theorem T, =
UlfQO
.
3.2, we can always find a unitary matrix T,: such that Tp
[
with a{ upper triangular and 7rl = -&‘,uo - &$qvo. Multiplying the 5 x 5 principal submatrix to the left by Z2 @ T,“, we get a unitary block Hessenberg matrix which can be written with -7; = CY,HZ_L~ - #y, v. and the g{ from above using block Schur parameters as -y;
5;
4 ?:
I To get a unitarily
12
1 Ii
-$ PF1 P,
Gr 1
similar matrix, we have to multiply the unitary block Hessenberg
matrix
M. Van Barel, A. Bultheel/Journal of Computational and Applied Mathematics 50 (1994) 545-563
552
to the right by I2 CDT, 63 Zm_4. Hence, we can repeat the same reasoning for Tl, T2, T3, . . . as for To. Before the final step of the algorithm is executed, we have the following unitarily similar (m + 1) x (m + 1) matrix: 12(m’-2)
12(r?*‘-2) -afit,-
G;G; . . .G’,,_,
-rum
P,“L,
8;tL2
-2
Vml-2
t&2
cd
-2
j&-2
’ LS I
A5
with S = 1 if m is odd, and 8 = 2 if m is even. The unitary matrices GL are built up using the prime block Schur parameters. The final step transforms the last (3 + 8) x (3 + 6) block
-+$-2 &‘-2
B =
$1:
,1
[’
;?:’
;;J
[7:;’
z:
z8].
[ l
If m is odd, there exists a unitary 2 x 2 matrix T,,r_l such that
with q,,,_, = -&_2,um~_2
- i;trf,_2ymj--l v,,,/-2 and al,,__, upper triangular. Finally we get
The last unitary matrix Gk,_, is built up by these block Schur parameters. l If m is even, there exists a 3 x 3 unitary matrix T,r_, such that TEL
[ g2;;,_2]
= [+I
)
with r,,,> _, as above and c$,, _, upper triangular. Finally we get
[
z2T;,_,] B [”
Tm,_l] = [ iFA’
These block Schur parameters
4. The solution
$1:
l]
[ z2 ;t
;v]
1
give us the final unitary matrices GL,_, and GL?,.
of the rational approximation
problem
After applying m steps of the algorithm constructed [u‘ OIT with H unitary block Hessenberg
above, we get Q”ZQ = H and Q” [f’
given by its block Schur parameters
- e’] =
yk, (Tk, &_kand Tk,
M. Van Barel, A. Bultheel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
k= 1,2,. . . , m’ - 1. For the moment we also assume that all gk are nonsingular, the regular case. If some flk is singular, we call it the singular case.
553
which we shall call
4.1. The regular case When the block Schur parameters gk are regular, the block Schur parameters allow us to define the orthonormal block polynomials +k in a recursive way. First, we need the following properties for the block columns of Q based on the same Schur parameters. Theorem 4.1. The block columns of Q satisfy the recursion ti!k+lck
=
ZQk
with initialization
+
Q;+&
&7’k,
Ql = Qi = [f’
= zQk7’: + Q;,
k=
1,2,3 ,...,
m’-
1,
- e’] a;‘.
Proof. From Q” [f’ - e’] = [(T: OIT, we get that Q, = [f’ - e’] ~0’. We set Q{ = Qt. From Q”ZQ = H, we get ZQ = QH = QGIG~.. . G,,! _, . Denote the (k + 1) st block column QGl G2Gk by Q:,, . Hence, by comparing the block columns of ZQ = QH, we get the recursion
of
Q;,, = Q:h + Qk+l%.
ZQk = -Q:?‘k + Qk+tak,
We get Qk+l out of the first equation and use this result to eliminate Qk+r from the second equation. Hence,
Q:,, = ZQkg;‘?k + Q;@-, +
Qk+l = (ZQk + Q;Yk)6’,
Ykg;‘?kh
Because [;T
;j
[;$
we can rewrite ~k’fk
$]=
[‘a I”i],
as $&L”
and &k + rkg;‘&
as 3,“.
This proves the theorem.
We define the block polynomials &, k = 0, 1,2, . . . , by a similar recursion block polynomials are orthonormal with respect to the scalar product (A(z),
B(z)),, =
•!
and show that these
~~TsB(z,,. i=l
Definition 4.2. We define the block polynomials based on the block Schur parameters: 4kcZ)u.k
= i$k-I
(z>
+
&&)Yk,
~~(z)~~=Z~k-I(Z)Y~+~:-,(Z)r
with the initialization Property 4.3.
C$kand c&, k = 0, 1,2, . . . , by the following recursion
k= 1,2 ,...,
m’-
(4.1)
1,
& = &, = u;‘.
The block polynomials
+k and 4; satisfy the following
properties
(with &,, = a:).
554
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
50 (1994) 545-563
l The polynomial qf~k has degree k, highest-degree coejficient u;‘u;’ ’ ’ ’ ~;?,a; and constant coefficient &;“LY;” +. . &;T, &F”@. a The polynomial q5: has degree k, highest-degree coejficient a;‘~;’ e. . a;‘Tk and constant coefj’icient a;“&,” . . . &iT, if?;“. l The block polynomials & and 4; are connected to the block columns Qk and Qi as follows:
Q;,, = D@s:,
Qk+l = D@kr with D = diag{ [fi jined similarly).
- e{]
, . . . , [f:,
-e;,]}
and
@k
=
[+k(zdT
4k(zdT
. . . +k(zrdTIT(@;
is
Proof. The properties immediately follow from the recurrence relations for the block polynomials and 4: and the fact that Yk, (7k, &_k,Tk are block Schur parameters. 0 The following result allows us to solve the discrete linearized problem on the unit circle. Theorem 4.4 (Orthonormality). The block polynomials respect to the scalar product (., *)u, i.e., (4i(z), @j(z)),
least-squares
&
rational approximation
C$kas dejined above are orthonormal =
de-
with
SijZz.
Proof. Using the results of the previous theorem, the scalar product can also be written as (41 (Z) 9rb;(Z >)c = @iDHD@; = Q,“,,Qj+I = SijZIv because QHQ = I,,,.
0
Let us now return to the rational approximation problem. To solve this problem, we shall write the solution in terms of the orthonormal block polynomials 4j. Because the highest-degree coefficients of these block polynomials are upper triangular, we can always write our solution n( z ) /d (z ) of degree LYas follows:
with (p,, 7,) # (0,O). Remember that the solution we are looking for has to be normalized. Property 4.3, we know that the highest-degree coefficient of 4a is equal to
From
(4.2) The specific value of x is not important for the moment. To normalize n( z ) /d( z ), we have to choose 7, = n,“, S2.j if n, # 0 or per = nj”, Sl,j if n, = 0. To solve the rational approximation problem, we have to minimize
M. Van Barel, A. BultheeUJournal
of Computational
and Applied Mathematics
SO (1994) 545-563
555
As we have seen in Section 2, we can make this entity smaller by putting &. and r/, equal to zero for (n(z),d(z)) andthedatato (lp,]*+]r,]*)‘/*. k=O, 1,. ..,a-l.Thisreducesthedistancebetween The way in which we have stated the rational approximation problem requires n(z ) /d( z ) to have degree cy. From the normalization conditions, it is clear that we should therefore take the “right” structure, i.e., 7, = I-I S2.i and pa = 0,
if IJ
and the “left” structure,
Pa =fi
Is2,jl6 I-J ISl.jI9 j=O
j=O
j=O
i.e.,
Si,j and 7, = 0,
if fi
IS*,jl3 fi
ISl,jl,
j=O
.i=O
.i=o
where, in case of equality of the two distances, we can choose between the two solutions. Note that this choice corresponds to our problem setting, but that we have the possibility to solve also slightly different problems in an equally simple way. For instance, we could have required n( z ) /d( z ) to be strictly proper, i.e., with IZ, = 0 and d, = 1. Then we had to choose the second solution, i.e., the one with 7, = 0. It is easy to generalize this to more complicated structures than the strictly proper one. Suppose for example that we look for a solution with ~1, = 0, IZ,_~ = 0 and d, = 1. From n, = 0, we easily derive that 7, = 0. From d, = 1, it follows that pa = nj”, sl,j. To make n,_i equal to zero, we have to take T,_~ equal to a specific value, possibly different from zero. We shall not follow this path in the sequel but concentrate our attention on the original problem setting. The other cases are quite similar. Let us investigate how good the rational approximation is. in terms of the parameters pk, ?-k, Theorem 4.5. If we write our rational approximant n(z)/d(z) k = 0,1,2 ,..., cy, then the approximation error made in the ith point z;, i = 1,2, . . . , m, can be measured as follows:
f,‘d(zi> - e(n(zi>= ~QG+I
[
:I] .
!A
Proof. This result follows immediately [j-i’ - e;]
@kc&)
=
QG+I.
the ith row of the block column Qk+i. For eXampk, if we take the Pn = n,“, Si,j. Th e approximation f:d(zi)
from the fact that
- eInCzi> = Qi,a+l
0
proper SOhtiOU, all pk and rk VaheS error in each point zi can thus be written as
Strktly
[
PU 0
I ’
i.e., the error is given by a scaled version of the first column of Qn+, .
are
zero
except
one,
M. Van Barel, A. Bultheel/Journal
556
of Computational and Applied Mathematics 50 (1994) 545-563
4.2. The singular case Until now, we have taken for granted that all the subdiagonal blocks uj were nonsingular. What will happen if one of the upper triangular matrices aj is singular, i.e., if si,j or s2,j is equal to zero?
Theorem 4.6. We shall assume that ga is the first block Schur parameter s2,a equal to zero. We de&e the block polynomial & as follows: #a = IfSl,a
=
Z&l
(2
> +
4&_,
(z
with ~1,~ or (if it exists)
ha.
0, then
gives us a normalized rational interpolant Otherwise, ~2,~ = 0 and
n( z ) /d( z ) with d, = 1.
gives us a normalized
n( z ) /d (z ) with n, = 1.
rational interpolant
Proof. Because the previous block Schur parameters compute
c$~ based on $a-1
and &_, . Moreover,
aj, j = 0, 1,2, . . . , a - 1, are regular, we can c$~ has degree (Y with the same highest-degree
coefficient as c$~_, , i.e., n,“;;’ fl,:‘. Using the relationship between the block columns of the unitary matrix Q and the orthonormal block polynomials @jy j = 0, 1, . . . , a - 1, we get
Qi,a+l u,=
[fi
-ejfja(Zi),
i=1,2
,...,
m,
with CT, =
[
Sl,o
S3,a
0
S2.a
If we summarize
1 *
the two cases as
we get that [fi
- ei]
tlz,’ [ I 1 =[fi’
- 414aC-G)
f
[I
=Qi,a+lua [f]
=O,
i=1,2
,...,
m,
for the specific choices of jj and ?. Hence, n( z)/d( z ) is a rational interpolant for the given data. By looking at the highest-degree coefficient, we see that this rational interpolant is normalized. This proves the theorem. 0
M. Van Barel, A. BultheeUJournal
of Computational
and Applied
Mathematics
50 (1994)
545-563
551
The rational function n( z ) /d( z ) from the previous theorem approximates the given data exactly, i.e., n( z)/d( z) is an interpolant of the given data. It is clear that this situation will appear rarely in practice when the data are only given with a certain finite precision.
5. Deriving a rational approximant
with real coefficients
In a lot of practical applications, e.g., digital filter design, model reduction, modal analysis, system identification, one looks for a rational approximant n( z ) /d( z ) with real coefficients. The algorithm, as described before, results in complex coefficients in general. If we want n( z)/d( z ) to have real coefficients, the data should be consistent with this restriction, i.e., if z is a given abscissa with corresponding estimated function value f/e, then f/a can be considered as an estimated function value in F. Using this fact, the initial data matrix M can be transformed using a simple unitary matrix into a similar real matrix. From this point on, all computations can be done using real arithmetic, resulting in real block Schur parameters, real block orthogonal polynomials and finally a real rational approximant. The simple unitary matrix to do this initial similarity transformation can be chosen for each data pair in several ways, e.g.,
because
6. Parallel implementation Regarding the amount of computational work, we can make a distinction between the computation of the unitary block Hessenberg matrix QHZQ from the data and the computation of the orthonormal block polynomials. Fixing the degree LY< m of the rational approximant implies that we can obtain the rational approximant from +a (z ) . Hence, it is not necessary to compute the complete block Schur parametrization but only the first part consisting of Gr , G2, . . . , G,. Compared to the complexity rn2 of computing the complete block Schur parametrization, this decreases the amount of computational work to ICY. The second part of the solution method requires the computation of &,+,, . . . , $J,, which requires a number of floating-point operations of the order of (Ye. Therefore, since (Y < m, it is important to decrease the execution time of the first part by parallel implementation. This could be done by a sort of pipeline method. Let us assume that we have cv < m processors. Each processor is coupled to a Gk, k = 1,2, . . . , a, matrix and updates the corresponding block Schur parameters in each step. Once the kth processor
M. Van Barel, A. Bultheel/Journal of Computational and Applied Mathematics 50 (1994) 54.5-563
558
IO-4
10-7 -
IO-10
IO-13IO-16
0
+ 10
5
1
IS
20
25
30
35
40
Fig. 1. The modulus of the function values of Example 7.1.
has computed CQ, Pk, &, & and Tk, the (k + 1)th processor can use this information to update the corresponding block Schur parameters. Once the pipeline is started up, the execution time is equivalent to O(m) floating-point operations.
7. Examples In this section, we give three examples. The first one is rather theoretical, the others are more practical. In the three examples, we look for a rational approximant with real coefficients (see Section 5). To measure the goodness of the rational approximant, we consider a relative and absolute approximation error in each point zi. We take all the weights wi(p) = 1. The absolute error that is plotted in the figures is then nothing else but Jrip)j. Similarly, when the relative error is plotted, we plot jr~P’/(fi/ej)l. The data is normalized in such a way that Cyl, (Ifil’ + lei12) = jlaoll; = 1 (11. IIF stands for the Frobenius norm). We can do this by computing Wi such that wi(lfi12
+
lei12)= $.
(7.1)
The value of the rational approximant we could also use the recurrence (4.1) the unit circle, as in the examples here, precision on a DEC workstation using
n( zi) /d( zi) is computed using the Homer scheme. Of course, where z is replaced by zi. For equally distributed points zi on one could also use FIT. All computations were done in double PRO-MATLAB with machine epsilon eps = 2.2204. 10-16.
Table 1
0
1 2 3 4 5
7.4349607 7.265 425 3 6.1803399. 6.1803399. 6.083 663 2 5.8550486.
’ 10-l . 10-l 10-l 10-l . 10-l 10-l
6.687403 1 5.5589297. 5.257311 1 6.305 926 4 5.473 828 7 . 5.123981 1.
10-l 10-l 10-l lo-l6 lo-l6 lo-l6
M. Van Barel, A. Bultheel/Journal
of Computational
and Applied Mathematics
50 (1994) 545-563
559
Table 2 Degree
Numerator coefficients
Denominator
0
3.152963 175 230421 IO-‘” 1.000000000000001 . loo 8.088600263652522. 10-l’ 1.000 000 000 000 000 10”
1.ooo 000 000 000 000 loo 1.633 234097 724558 . lo--l6 1.576481587615210. lo-l6 1.435 322 868 745 093 1O-l6
1 2 3
coefficients
Example 7.1. In the first example, we consider the m = 40 interpolation m, with function values i= 1,2,..., fi
= Z” + ’
q
Z.
points zi = exp(ij27r/m),
”
i.e., the polynomial z (Z - j) (z + j) evaluated at the points zi. The modulus of these function values is plotted in Fig. 1. The errors connected to the “left” and “right” solution of the linearized problem are given in Table 1. From this table, it is clear that we should choose cx equal to 3 and take the “right” structure for the rational approximant. Table 2 shows the coefficients of the numerator and denominator polynomial of the rational approximant of degree LY= 3. The absolute and relative errors are indicated in Fig. 2. We can say that the rational approximant we obtained here is very good, which, for this example, of course will not come as a surprise. In the next example circle.
we try to find a good approximant
for the exponential
function
on the unit
Example 7.2. The points zi are the same as in the first example. The function values are fi/ei = eXp(Zi), i = 1,2,. , ., m. When looking at the errors of the solutions of the linearized problem, we choose to take a “right” solution of degree 3. The absolute and relative errors are given in Fig. 3. By executing one iteration step of the Loeb-Wittmeyer algorithm adapting the weights w, as described in Section 2, we get the errors of Fig. 4. Note that, using the normalization of (7.1)) the relative error is almost a constant over the set of points zi. However, doing one iteration step, it is the absolute error having this property. A similar phenomenon can also be observed in the next example. Moreover, the winding number of the absolute error function exp( z ) - n( z ) /d( z ) on the unit circle
IO--
* 0
lo-)2 5
IO
15
20
25
30
35
40
* 0
5
10
15
Fig. 2. The absolute and relative errors of Example 7.1.
2”
25
30
35
1
40
M. Van Barel, A. BultheeUJournal of Computational and Applied Mathematics 50 (1994) 545-563
560
Fig. 3. The absolute and relative errors of Example 7.2.
is equal to deg IZ+ deg d + 1. In several examples, the error function for a best Chebyshev rational approximant n//d’ closely approximates a perfect circle about the origin having a winding number deg n’ + deg d’ + 1 [ 131. Hence, the rational approximant obtained here after doing one iteration step will be a close candidate to a best Chebyshev rational approximant on the unit circle. A more practical example is the following one. Example 7.3. As in the previous examples, we assume that the points zi are equally distributed on the unit circle, i.e., we take zi = exp(ij27r/m), i = 1,2, . . . , m, with m = 200. We choose in a certain random way five zeros with modulus between 0 and 0.99. Similarly, we choose three poles with modulus between 0.90 and 0.99 and two poles with modulus between 0 and 0.90. With these poles and zeros, a rational function is constructed with real coefficients. Hence, we add the complex conjugates of the poles and the zeros. The data fi/ei is equal to this rational function evaluated in zi, but we have added a uniformly distributed absolute error of 10P4. The modulus of the function values is shown in Fig. 5. Examining the errors connected to the linearized problem, we choose for a “right” solution of a reduced degree 8. The absolute and relative errors are given in Fig. 6. Doing one iteration as described in Section 2, we obtain the absolute and relative errors given in Fig. 7. The
+++*+*+*+**++*++++++++**+*+++**+*-
IO-J _f.ff.f
+
IO-J
*
* f
.
f
*+ +**
IO.4 0
5
10
I5
20
2s
30
35
40
Li7rT-J 0
Fig. 4. The absolute and relative errors of Example 7.2.
+
*
M. Van Barel, A. Bultheel/Journal
of Computational and Applied Mathematics 50 (1994) 545-563
561
102
++
: +
t
t 10’
* +++
5 * *+
100
/
* ; *
i
*
t
+ t +* + t* f* *f
t ** f : c
w++++
Fig. 5. The modulus of the function values of Example 7.3.
left side of Fig. 8 shows the poles (indicated by *) and the zeros (indicated by o) of the original rational function of degree 10 in the complex plane. The right side shows the poles and zeros of the rational approximant after doing one iteration step. Note that the dominant poles and zeros, i.e., the ones closely to the unit circle, are well approximated. These three examples show that the rational approximant approximates the given data quite well in the proper rational sense. Moreover, applying only one iteration step of the algorithm of LoebWittmeyer already gives an absolute error whose magnitude does not change very much anymore over the set of points zi. Because the amount of computational work is small (O(ma) and, in parallel, O(m) floating-point operations), the designer of a rational approximant has a lot of freedom in choosing the degree, the iteration method (how are the weights adapted), the number of iteration steps, . . . . The errors connected to the “left” and “right” solutions of the linearized problem give a good indication for the designer to choose between the several possible structures for the rational approximant (degree, “left” or “right”) . In the context of discrete-time linear systems, we can interpret the rational approximation problem solved in this paper as a system identification problem. The data fi/ei can be considered as
I
IO-5 ’
0
20
40
60
80
100
120
140
160
180
200
0
20
I 40
60
x0
Fig. 6. The absolute and relative errors of Example 7.3.
100
I20
140
I60
180
200
M. Van Barel, A. BultheeUJournal of Computational and Applied Mathematics 50 (1994) 54.5-563
562
Fig. 7. The absolute and relative errors of Example 7.3.
measurements of the frequency response of a discrete-time system. We are now looking for a linear model (a rational function) of limited complexity (small degree) such that the measured frequency response and the frequency response of the model almost coincide (taking into account the error on the measurements). For an introduction to system identification one could look at, e.g., [7, Chapter 81 and [ 11, Chapter 61. The rational approximation problem can also be seen as a model reduction problem. In the last example, we start with a model of degree 10 and we approximate this by a model of degree 8 [ 5, Section 191.
8. Conclusion We have developed an efficient method to solve the discrete linearized least-squares rational approximation problem. We have shown that this method can be implemented in a pipelined way on a parallel architecture. The numerical examples show that the computed rational approximant is also quite good as a solution of the proper rational approximation problem. It becomes a very good proper rational approximant after only one iteration step of the algorithm of Loeb-Wittmeyer. Because the algorithm developed here is so efficient, it could be used as the kernel of a graphical
-1
-0.5
0
0.5
I
Fig. 8. The poles and zeros of the original and the approximating
rational function of Example 7.3.
M. Van Barel, A. Bultheel/Journal
of Computational
and Applied Mathematics
50 (1994) 54.5-563
563
interactive program to design rational approximants for large data sets. As domains of application, we mention digital filter design, modal analysis, model reduction, system identification, . . . . In future work, this algorithm will be generalized to other structures for the rational approximant, i.e., the maximum degree of numerator and denominator polynomial can be chosen freely. Also more than two data columns f’ and e’ will be allowed, leading to block orthonormal polynomials with the dimension of the blocks bigger than 2.
References ‘[l] G.S. Ammar and W.B. Gragg, O(n*) reduction algorithms for the construction of a band matrix from spectral data, SIAM _I. Matrix Anal. Appl. 12 ( 199 1) 426-43 1. [ 21 G.S. Ammar, W.B. Gragg and L. Reichel, Constructing a unitary Hessenberg matrix from spectral data, in: G.H. Golub and P. Van Dooren, Eds., Proc. NATO Advanced Study Institute on Numerical Linear Algebra, Digital Signal Processing and Parallel Algorithms, Leuven, 1988 (Springer, Berlin, 199 1) 385-395. [ 31 1. Barrodale and J.C. Mason, Two simple algorithms for discrete rational approximation, Math. Comp. 24 ( 112) (1970) 877-891. [4] D.L. Boley and G.H. Golub, A survey of matrix inverse eigenvalue problems, Inverse Problems 3 (1987) 595-622. [5] A. B&heel and M. Van Barel, Padt techniques for model reduction in linear system theory: a survey, J. Comput. Appl. Math. 14 (3) (1986) 401-438. Updating and downdating of orthogonal polynomials with data fitting [6] S. Elhay, G.H. Golub and J. Kautsky, applications, SIAM J. Matrix Anal. Appl. 12 (1991) 327-353. [7] P. Faurre and M. Depeyrot, Elements of System Theory (North-Holland, Amsterdam, 1977). [8] W.B. Gragg and W.J. Harrod, The numerically stable reconstruction of Jacobi matrices from spectral data, Numer. Math. 44 ( 1984) 317-335. [ 91 L. Reichel, Fast QR decomposition of Vandermonde-like matrices and polynomial least squares approximation, SIAM .I. Matrix Anal. Appl. 12 ( 1991) 552-564. [lo] L. Reichel, G.S. Ammar and W.B. Gragg, Discrete least squares approximation by trigonometric polynomials, Math. Comp. 57 ( 1991) 273-289. [ 111 W.J. Rugh, Mathematical Description of Linear Systems (Marcel Dekker, New York, 1975). [ 121 H. Rutishauser, On Jacobi rotation patterns, in: N.C. Metropolis, Ed., ExperimentalArithmetic, High Speed Computing and Mathematics, Proc. Sympos. Appl. Math. 15 (Amer. Mathematical Sot., Providence, RI, 1963) 219-239. [ 131 L.N. Trefethen, Near-circularity of the error curve in complex Chebyshev approximation, _/. Approx. Theory 31 ( 1981) 344-367. [ 141 M. Van Bare1 and A. Bultheel, A new approach to the rational interpolation problem, J. Comput. Appl. Math. 32 (l&2) (1990) 281-289. [ 151 M. Van Bare1 and A. Bultheel, A new approach to the rational interpolation problem: the vector case, J. Comput. Appl. Math. 33 (3) ( 1990) 33 l-346. [ 161 M. Van Bare1 and A. Bultheel, A parallel algorithm for discrete least squares rational approximation, Numer. Math. 63 (1992) 99-121.