# Approximate Fekete points for weighted polynomial interpolation.

1. Introduction. In the framework of polynomial interpolation, Fekete points are points that maximize the Vandermonde determinant (in any polynomial basis) on a given compact set and thus ensure that the corresponding Lebesgue constant grows (at most) algebraically, being bounded by the dimension of the polynomial space. Their analytical properties and their efficient computation are still essentially open research problems in the multivariate setting; see, e.g., [11, 44] and references therein. In particular, the computation of Fekete points requires solving large scale optimization problems already at moderate degrees. Much more is known in the univariate case, but the computational problem is still open in one complex variable (where, however, good alternatives are known, like Fejer or Leja-like points; see [3, 27]).

In some recent papers [9, 10, 39], a greedy algorithm has been studied, that computes (multivariate) approximate Fekete points by extracting maximum volume submatrices from rectangular Vandermonde matrices on suitable discretization meshes. It works on arbitrary geometries and uses only optimized tools of numerical linear algebra (essentially QR-like factorizations). There is a strong connection with the theory of admissible meshes for multivariate polynomial approximation, recently developed by Calvi and Levenberg [15]. There are also good perspectives in the application to numerical cubature and to the numerical solution of PDEs by collocation and discrete least squares methods [45]. A renewed interest is indeed arising in methods based on global polynomial approximation; see, e.g., [29].

The algorithm can be described in a very general functional, not necessarily polynomial, setting. Given a compact set k [member of] [R.sup.d] (or [C.sup.d]), a finite-dimensional space of linearly independent continuous functions,

(1.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.],

and a finite set {[[xi].sub.1], ..., [[xi].sub.N]} [member of] K, we can construct the Vandermonde-like matrix,

V([[xi].sub.1], ..., [[xi].sub.N]) = [[v.sub.ij]] := [[phi].sub.j] ([[xi].sub.i])].

If det(V([[xi].sub.1], ..., [[xi].sub.N]) [not equal to] then the set {[[xi].sub.1], ..., [[xi].sub.N]} is unisolvent for interpolation in [S.sub.N]

and

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

is a cardinal basis, i.e., [[PSI].sub.k] ([[xi].sub.k]) = [[delta].sub.jk] and

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

interpolates any function f at {[[xi].sub.1], ..., [[xi].sub.n]}. In the case that such points maximize the (absolute value of the) denominator of (1.2) in [K.sup.N] (Fekete points), then [parallel] [[PSI].sub.j] [[parallel].sub.[infinity]] [less than or equal to] 1 for every j, and thus the norm of the interpolation operator [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is bounded by the dimension of the interpolation space,

(1.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Clearly, Fekete points, as well as the "Lebesgue constant" [[LAMBDA].sub.N], are independent of the choice of the basis in [S.sub.N], since the determinant of the Vandermonde-like matrices changes by a factor independent of the points (namely the determinant of the transformation matrix between the bases).

The maximization in [K.sup.N] is a nonlinear optimization problem in dN real (or complex) variables. The idea of the algorithm is to maximize on a suitable discretization mesh of K,

(1.5) X = {[x.sub.i]} [member of] K, 1 [less than or equal to] i [less than or equal to] M, M > N,

i.e., to construct the rectangular M x N Vandermonde-like matrix

V([x.sub.i], ..., [x.sub.m]) = [[v.sub.ij]] := [[[phi].sub/j]([x.sub.i])],

and to extract from it a maximum volume N x N square submatrix. Observe that such a discrete nonlinear optimization problem is known to be NP-Hard (cf. [17]), but an approximate solution can be obtained by the following greedy algorithm applied to A = [V.sup.t]:
```Algorithm greedy (max volume submatrix of A [member of]
[R.sup.NxM], M > N)

* ind= [] ;

* for k = 1, ..., N

-- "select the largest norm column [col.sub.ik] (A)"; ind =
[ind, [i.sub.k]];

-- "remove from every column of A its orthogonal projection
onto [col.sub.ik]";
end;
```

which works when A is full rank, and gives the set of indexes ind = (i 1,..., iN) corresponding in our problem to the approximate Fekete points

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

The algorithm can be conveniently implemented by the well-known QR factorization with column pivoting, originally proposed by Businger and Golub in 1965 [12], and used for example by the MATLAB "mldivide" or "\" operator in the solution of underdetermined linear systems (via the LAPACK routine DGEQP3; cf. [25,31]). The full algorithm proposed in [39], applied in the present general setting, can be summarized in a MATLAB-like notation as follows:
```Algorithm AFP (Approximate Fekete Points by iterative refinement)

* take a suitable discrete subset X = ([x.sub.1], ..., [x.sub.M])
[member of] K, M > N;

* [V.sub.0] = V([x.sub.1], ...,[x.sub.m]); [P.sub.0] = I;

* for k = 0, ..., s - 1

[V.sub.k]c = [Q.sub.k][R.sub.k]; [U.sub.k] = inv([R.sub.k]);
[V.sub.k+1] = [V.sub.k][U.sub.k]; [P.sub.k+1] =
[P.sub.k][U.sub.k];

end;

* A = [V.sup.t.sub.s]; b = [(1, ..., 1).sup.t]; (the choice of
b is irrelevant in practice)

* w = A\b; ind = find(w [not equal to] 0); [X.sub.*] = X(ind);
```

The greedy algorithm above is implemented directly by the last row of Algorithm AFP (in MATLAB), irrespectively of the actual value of the vector b. The for loop above implements a change of polynomial basis from ([[phi].sub.1], ..., [[phi].sub.N]) to the nearly-orthogonal basis (q.sub.1], ..., [q.sub.N]) = (([[phi].sub.1], ..., [[phi].sub.N]) [P.sub.s] with respect to the discrete inner product (f, g) = [SIGMA] f([x.sub.j])[bar.x]([x.sub.j])> whose purpose is to overcome possible numerical rank-deficiency and severe ill-conditioning arising with nonorthogonal bases (usually s = 1 or s = 2 iterations suffice); for a complete discussion of this algorithm we refer the reader to [ 9, 10, 39].

REMARK 1.1. The effectiveness of Algorithm AFP in producing good interpolation points depends on the distribution of the initial discretization points, which in turn should take into account the geometry of the domain as well as the peculiarity of the function space. In the case of total-degree nonweighted polynomial interpolation, [S.sup.N] = [P.sup.d.sub.n], N = dim([P.sup.d.sub.n](K)), it is known that good starting meshes are the (weakly) admissible meshes studied in [15], since approximate Fekete points extracted from such meshes have good interpolation properties and the same asymptotic behavior of the true Fekete points; cf. [9, 10, 39]. In particular, the associate discrete measure [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] converges weak-* to the equilibrium measure of K in the sense of (pluri)potential theory [9].

REMARK 1.2. Observe that if we take b = ([m.sub.1], ..., [m.sub.N]).sup.t][P.sub.s] in the Algorithm AFP, where [m.sub.j] = [[integral].sub.K] [[phi].sub.j]- (x)d/[mu] are the "moments" of the original basis with respect a given measure on [mu], then w(ind) = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is an array of weights of a cubature formula which is exact on [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. For example, in the total-degree polynomial case, [S.sup.N] = [P.sup.d.sub.n](K), with d[mu] = dx, the moments of any polynomial basis can be computed over arbitrary geometries by the Gauss-Green cubature formula, based on spline tracking of the boundaries, developed in [40].

2. Weighted polynomial interpolation. The literature on weighted polynomial interpolation is very extensive, since it concerns a variety of theoretical and applied topics, such as rational interpolation with prescribed poles, weighted potential theory, Gaussian quadrature, numerical treatment of integral equations, design of weighted digital filters in signal processing, and many others. Extensive references to such a vast literature are beyond the scope of this work; we suggest that the interested reader consult, for example, [22, 26, 30, 36] with the references therein.

In the present paper, which is mainly of a computational and experimental character, we focus on approximating Fekete points for weighted polynomial interpolation on compact sets, in the following framework; cf. (1.1)-(1.3). Let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where k [member of] [R.sup.d] (or [C.sup.d]) is a compact set and ([P.sub.j]) is a basis of the total-degree polynomial space,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We have in mind two distinct situations

* prescribed singularities: we interpolate in [S.sub.N] functions of the form

f = wg

where g is regular and the weight function absorbs the singularities of f, like, e.g., real or complex poles not belonging but possibly close to K.

* weighted norms: we interpolate directly a regular function g by a polynomial p [member of] [P.sup.d.sub.n](K), but the error is measured with a weighted norm,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where weights in different ways different parts of the domain, e.g., the union of two disjoint intervals, which is relevant to the design of digital filters.

As will be shown in the examples, the important thing with prescribed singularities is to adjust the basis according to the weight function. Once a suitable weighted basis for interpolation is chosen, weighted approximate Fekete points give slightly better results than nonweighted ones.

On the other hand, the two instances above are two faces of the same coin, as is shown by the following observation. Let the set {[[xi].sub.1], ... [[xi].sub.N]} be unisolvent for interpolation in [S.sup.N]. Then, it is also unisovent for interpolation in [P.sup.d.sub.n](K), in view of the following relation between Vandermonde-like matrices,

(2.1) [w([[xi].sub.i])[p.sub.j]([[xi].sub.i])] = diag(w([[xi].sub.i])) [[p.sub.j](&)].

Moreover, it is easy to see that the cardinal functions for interpolation in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] where [l.sub.j] is the fundamental Lagrange polynomial for the point [[xi].sub.j]. Observe that unisolvence implies that w([[xi].sub.j]) [not equal to] 0 for all j. On the other hand, the existence of unisolvent interpolation sets for [S.sub.N], for example weighted Fekete points, is guaranteed by continuity of w as soon as supp(w) [intersection] K is polynomial determining (i.e., any polynomial vanishing there is identically zero).

Concerning the interpolation operators, it follows immediately by uniqueness that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

and thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Concerning convergence, it is worth observing that, at least with "exact" Fekete points for [S.sub.N],

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

i.e., convergence is certainly guaranteed as soon as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

(we stress, however, that often the Lebesgue constant of (approximate) Fekete points has a much slower growth than that estimated by (1.4)). In this regard, Jackson-like theorems on the rate of convergence of best polynomial approximation on K are needed; cf., e.g., [33]. In one complex variable, for example, by a classical result of Walsh and Russell [43], we have convergence for every function g analytic in a compact set K with connected complement, even when K = [[union].sup.m.sub.i=1] [K.sub.i] is a finite disjoint union of compacts with Jordan boundary, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], are possibly different analytic functions. Observe that if w is analytic itself in K, the maximum modulus of the determinant of the weighted Vandermonde matrix (2.1) in [K.sup.N], is attained on the product of boundaries [([partial derivative]K).sup.N] [subset of equal to] [partial derivative][K.sup.N] by the maximum principle applied to each variable. Indeed, it can even be proved that the Shilov boundary of [K.sup.N], i.e., the smallest subset of the topological boundary d[K.sup.N] where every holomorphic function of N complex variables attains its maximum modulus, is contained in [([partial derivative]K).sup.N]; cf., e.g., [37, [section]2.5]. Thus we can compute the Fekete points on the boundary [partial derivative]K of K. This is an advantage, since also geometrically we deal with a one-dimensional instead of a two-dimensional problem.

Before presenting the numerical examples, it is worth discussing briefly the following problem: What is a reasonable distribution of the starting mesh X, from which we can extract approximate Fekete points by algorithm AFP? In the case of nonweighted polynomial interpolation, a guideline is given by the theory of "admissible meshes" for polynomial approximation recently developed by Calvi and Levenberg in [15]. The key feature of an admissible mesh X is a polynomial inequality like

(2.3) [[parallel]p[parallel].sub.k] [less than or equal to] C [parallel]p[parallel]x, p [member of] [P.sup.d.sub.n](K),

which ensures that the Lebesgue constant of Fekete points of X can be bounded proportionally to that of exact Fekete points of K,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

since for the fundamental Lagrange polynomials we have the bound [parallel][l.sub.K] [less than or equal to] C [parallel] [l.sub.j] [[parallel].sub.x] [less than or equal to] < C; see [15, [section]4.4]. Observe that necessarily card(X) [greater than or equal to] dim([P.sup.d.sub.n](K)).

Starting from (2.3), a rough functional inequality can be obtained also in the space [S.sub.N] = w[P.sup.d.sub.n](K), at least if w [not equal to] 0 in K. Indeed, if X is an admissible mesh for nonweighted polynomial interpolation,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

for suitable [??] [member of] K, [eta] [member of] X. Thus we get, if w [not equal to] 0 in K,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

i.e., X is an admissible mesh also for weighted polynomial interpolation. An estimate on the growth of the Lebesgue constant of Fekete points of X in the space [S.sub.N] now follows,

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

We observe, however, that (2.4) turns out to be largely an overestimate of the actual growth of the Lebesgue constant. In order to get more suitable meshes and more refined bounds one should take into account the specific structure of the weight function w, possibly by using Markov-like inequalities for weighted polynomials, to mimick the construction of [15, Thm. 5].

A possible refinement can be given when [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (the discrete subsets [X.sub.i] are not disjoint, in general). In this case we obtain easily the estimate,

(2.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

3. Numerical examples. The study of rational interpolation and approximation of univariate analytic functions with prescribed poles has a long history, dating back to the fundamental work by Walsh and successive extensions by Bagby; cf. [2,42]. Also its computational issues have been deeply investigated; see, e.g., [4, 5, 16, 34, 41] and references therein.

Here we begin by considering weighted polynomial interpolation in the space [S.sub.N] = w[P.sub.n]([-1,1]) with a weight function

w(x) = 1/[[pi].sub.m](x), [[pi].sub.m] [member of] [P.sub.m]([-1,1]);

where the real and complex zeros of [[pi].sub.n] do not belong (but possibly are close) to the interval.

In the literature computational methods have been studied to find almost optimal points for rational interpolation with prescribed poles. For example, one can compute and use as interpolation points in [S.sub.N] the zeros of the monic polynomial [p.sub.n] [member of] [P.sub.n] such that the ratio [p.sub.n]/[[pi].sub.m] has the smallest possible max-norm (a min-max approach like that leading to Chebyshev polynomials and Chebyshev points for polynomial interpolation). This can be done by transforming the problem to a numerical eigenvalue problem; cf., e.g., [41] and references therein.

The first four examples below show that the computation of approximate Fekete points for [S.sub.N] is a natural and effective alternative approach to the problem of rational interpolation with prescribed poles in one real variable. On the other hand, the algorithm is general purpose, and is able to handle even complex or bivariate instances, as well as other kinds of singularities. It can also be applied to the construction of weighted real or complex polynomial filters useful in signal processing.

EXAMPLE 3.1. (real interval, one real pole) We interpolate the function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

which has a real pole of the second order at x = 1 + [epsilon].

In Figure 3.1 we compare polynomial interpolation at the Chebyshev-Lobatto points, with interpolation in w[P.sub.n] again at the Chebyshev-Lobatto points for

w(x) = [w.sub.[epsilon]](x) = [(1 + [epsilon] - x).sup.-2],

and at approximate Fekete points computed by Algorithm AFP using the Chebyshev polynomial basis (i.e., [[phi].sub.j] (x) = [w.sub.[epsilon]](x)[T.sub.j](x), [T.sub.j](x) = cos (j arccos x), j = 0, ..., n), with 2 refinement iterations and a starting mesh X of 1000 equispaced points; cf. (1.5). The density of X ensures that it is an admissible mesh in the sense of (2.3) for (nonweighted) polynomial interpolation up to the highest degree (n = 30), since it has stepsize h < 2/[n.sup.2]. This result relies on the classical Markov polynomial inequality [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]; cf. [6] and the construction in [15, Thm. 5].

In this case, estimate (2.4) has no practical usefulness, since it becomes [LAMBDA]n = O([[epsilon].sup.-2])N = O(n [[epsilon].sup.-2]), which is several orders of magnitude above the actual growth.

[FIGURE 3.1 OMITTED]

To improve the construction, one should have at hand a tight Markov-like inequality for rational functions of the form P/[[pi].sub.m], p [member of] [P.sub.n] and m fixed, whereas this subject seems to have been studied in the literature for the case of ( n, n) rational functions (with denominator of degree exactly n).

Some improvement in the estimate of the Lebesgue constant, however, can obtained following (2.5). Given the mesh X with constant stepsize as above, let us take [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] Then K = [union] [K.sub.i], X = [union] = [X.sub.i], and all the [X.sub.i] are admissible meshes with constants [C.sub.i] [less than or equal to] C. On the other hand, it is not difficult to show that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] Since h = O([n.sup.-2]), the mesh X being an admissible mesh for nonweighted polynomial interpolation, we get [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and by (2.4)-(2.5), being N = n + 1, the final estimate [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.].

EXAMPLE 3.2. (real interval, two conjugate poles) We consider here the function

f(x) = cos (x)/[[epsilon].sup.2] + [x.sup.2], x [member of] K = [-1, 1],

which has two complex conjugate poles of the first order at x = [+ or -]ie, taking

[w.sub.[epsilon]] (x) = [([[epsilon].sup.2] + [x.sup.2]).sup.-1].

The numerical results are collected in Figure 3.2. Observe that, in both Examples 3. 1 and 3.2, the choice of the interpolation space is much more relevant than the choice of the points.

[FIGURE 3.2 OMITTED]

[FIGURE 3.3 OMITTED]

Nevertheless, weighted polynomial interpolation at approximate Fekete points shows small Lebesgue constants, close to those of nonweighted polynomial interpolation (growth like O(logn)), and the best interpolation errors also for small [epsilon]. It is interesting to have a look at the distribution of the interpolation points; see Figures 3.3 and 3.4. In Figure 3.4 we show the distribution function of the interpolation points, that is the fraction of points which are smaller than or equal to x. Notice that for [epsilon] = 1 the approximate Fekete points distribute like the Chebyshev-Lobatto points, whereas for smaller they also tend to cluster at the point in the interval nearest to the pole(s).

EXAMPLE 3.3. (two disjoint intervals, one real pole in between) In the case of a single real interval, other algorithms exist to compute near-optimal points for rational interpolation with prescribed poles, like that described in [41]. On the other hand, one of the strengths of Algorithm AFP is that it can work on quite general compact sets, with real or complex variables. We consider now the case of simultaneous interpolation on two disjoint real intervals,

f(x) = cos(x)/[(x + [epsilon]).sup.2], [member of] K [-1,-0.6][union][0,1],

[FIGURE 3.4 OMITTED]

where f has a second order pole at x = -[epsilon], which lies between the two intervals for 0 < [epsilon] < 0.6. We take

[w.sub.[epsilon]](x) = [(x + [epsilon]).sup.-2],

and we discretize the first interval by 200 and the second by 500 equispaced points, respectively. This gives an admissible mesh for (nonweighted) polynomial interpolation on K up to degree n = 30, as union of two admissible meshes; cf. [15] and the discussion in EXAMPLE 3.1. The interpolation spaces are [P.sub.n] with approximate Fekete points for nonweighted polynomial interpolation, wePn again with the latter points, and [w.sub.[epsilon]][P.sub.n] with approximate Fekete points for weighted polynomial interpolation (basis [w.sub.[epsilon]][T.sub.j]} in Algorithm AFP). Figure 3.5 collects the numerical results, that are comparable to those of EXAMPLE 3.2.

In Figures 3.6 and 3.7 we show the distribution of the interpolation points. In the case of nonweighted interpolation, the physical interpretation is that we have computed an approximate equilibrium configuration of N repelling equal charges located on the two intervals. We observe that there are more interpolation points in the largest interval, and that they tend to cluster at the intervals endpoints, more rapidly at the external ones and in case of weighted interpolation also at the endpoint nearest to the pole. This has been confirmed in several other numerical experiments (by increasing n and decreasing [epsilon]), which are not reported for brevity's sake.

EXAMPLE 3.4. (complex disk, one real pole) In order to show the flexibility of the method, we consider again the function of EXAMPLE 3.1, but now on the complex unit disk

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

[FIGURE 3.5 OMITTED]

As already observed, since the weight function [w.sub.[epsilon]] (x) = [(1 + [epsilon] - x).sup.-2] is analytic in a neighborhood of K, we can compute the Fekete points on the boundary [partial derivative]K = {x [member of] C : |x| = 1} by the maximum principle. Fekete points for nonweighted polynomial interpolation are one of the few explicitly known cases, namely any sequence of N = n +1 equispaced points on the unit circle is a set of Fekete points for interpolation degree n. The approximate Fekete points for weighted polynomial interpolation have been computed by Algorithm AFP using the monomial basis and two refinement iterations, starting form a mesh of 1000 equispaced points on the unit circle.

The numerical results are reported in Figure 3.8, and in Figure 3.9 the approximate Fekete points for degree n =15, corresponding to [epsilon] =1 and [epsilon] =0.01, are displayed. Notice that for the smallest value of [epsilon] the points tend to cluster at x = 1, the point of the circle nearest to the pole x = 1 + [epsilon].

EXAMPLE 3.5. (two real variables in the square, line of algebraic singularities) We give now an example in two real variables, taking the function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

This function is analytic up to an entire line of singularities at x(1) = 1 + e.

[FIGURE 3.6 OMITTED]

[FIGURE 3.7 OMITTED]

[FIGURE 3.8 OMITTED]

[FIGURE 3.9 OMITTED]

In Figure 3.10 we show the interpolation errors and the Lebesgue constants corresponding to different interpolation spaces and nodes, and in Figure 3.11 we plot the approximate Fekete points corresponding to degree n = 10. As a starting mesh to extract approximate Fekete points we have taken a 120 x 120 uniform grid, which is an admissible mesh for (nonweighted) polynomial interpolation up to n =10 (being the product of two one-dimensional admissible meshes; cf. [15]). The comparison is with nonweighted interpolation at the so-called "Padua points", the first known example of nearly optimal points for total degree polynomial interpolation in two variables, with a Lebesgue constant increasing like log squared of the degree; cf. [7, 8, 14].

EXAMPLE 3.6. (two disjoint intervals, nonalgebraic singularities) Algorithm AFP is general purpose, and can handle without problems even nonalgebraic singularities, for example singularities of derivatives inside the domain. A nontrivial example of this kind is given by the following function, defined on the union of two disjoint intervals K = [-1, -0.6] [union] [0,1],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Taking as a weight function,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

the function f is factorizable as f = wg with g (separately) analytic in each of the two intervals (since the Maclaurin series of sin and arctan have only odd powers), whereas w absorbs the singularities of f (which is only Holder continuous). The numerical results are collected in Figures 3.12 and 3.13. Notice that the errors of weighted polynomial interpolation are larger than in EXAMPLE 3.4, since here we deal with different analytic functions on disjoint intervals, and the convergence rates of best polynomial approximations are different.

EXAMPLE 3.7. (real digital filters) Simply stated, in the design of weighted polynomial digital filters of FIR (Finite Impulse Response) type, one seeks a polynomial p of degree n that approximates a function g, termed the response of the filter, on a real or complex compact set K, in the sense that the weighted norm [parallel]w(p-g)[[parallel].sub.K] is small, where w is a suitable continuous and nonnegative weight function. In many applications K is a finite union of disjoint intervals or complex arcs (the frequency bands), and both g and w are piecewise constant functions (the values 0 and 1 of g corresponding to the so-called stopband and passband, respectively).

[FIGURE 3.10 OMITTED]

[FIGURE 3.11 OMITTED]

We cannot even try to give an appropriate quoting of the enormous literature on this important subject of signal processing. To give only some highlights, it is worth noting the classical paper [32], where the popular Parks-McClellan algorithm for the design of optimal real equiripple filters, based on the Remez exchange algorithm, was originally proposed, as well as the treatise [13]. The method has then been extended in various directions to deal also with complex filters; see, e.g., [26] and references therein. Other methods are based, for example, on potential theory and conformal mapping; see, e.g., [21, 23, 38]. The use of polynomial filters is also interesting within numerical linear algebra; see, e.g., [24, 35].

[FIGURE 3.12 OMITTED]

[FIGURE 3.13 OMITTED]

Since we use here for the first time algorithm AFP for the construction of polynomial filters based on interpolation at approximate Fekete points, we begin by a simple nonweighted case, an example of a high-pass filter. We consider the response function (this example is taken from [21, [section]6]),

g(x) = {0, x [member of] [K.sub.1] = [-1, -0.4], 1, x [member of] [K.sub.2] = [-0.3,1],

defined on K = [K.sub.1] [union] [K.sub.2], where [K.sub.1] is the stopband and [K.sub.2] is the passband (the interval (-0.4, -0.3) is the transition band in signal processing terminology). As weight function we take w = 1, and as initial discretization an admissible mesh obtained by union of two uniform admissible meshes of [K.sub.1] and [K.sub.2], respectively.

The numerical results are collected in Figure 3.14. The amplitude of the oscillations near the internal endpoints (a sort of Gibb's phenomenon well studied in [38]) depends on the length of the transition band. In Figure 3.15, for comparison, we show the filter of degree n = 30 with a transition band of double length, namely for [K.sub.2] = [-0.2,1], and the interpolation errors up to degree n = 60.

We observe that the quality of the filters is lower than that of optimal equiripple filters (but not as much as could be predicted by the Lebesgue constants, cf. (2.2)), and even slightly lower than that of filters obtained via numerical conformal mapping; cf. [21, [section]6]. This means that interpolation at approximate Fekete points cannot be considered a real competitor in standard instances of polynomial filtering. On the other hand, its strength consists in a higher flexibility, that makes it applicable to quite general compact sets K, response functions g and weight functions w, in one or even several real or complex variables.

In order to give an example of a weighted filter, we consider the following multiband response function on the union of three disjoint intervals K = [K.sub.1] [union] [K.sub.2] [union] [K.sub.3],

g(x){0, x [member of] [K.sub.1] = [-1, -0.4], 1, x [member of] [K.sub.2] = [-0.2, 0.2], 0.5, x [member of] [K.sub.3] = [0.4, 1],

with a piecewise constant weight function corresponding to a triple of positive weights ([w.sub.1], [w.sub.2], [w.sub.3])

w(x) = [w.sub.i], x [member of] [K.sub.i],i = 1, 2, 3.

This is only an illustrative example, but it is worth recalling that multiband polynomial filters are a standard tool in digital signal processing; cf. [13].

In Figure 3.16 we show four multiband filters, obtained by interpolation at approximate Fekete points of degree n = 30, corresponding to different choices of the weights. In Figure 3.17 we plot the estimated Lebesgue constants, up to degree n = 30. Again, the starting discretization mesh is obtained by taking the union of three admissible meshes for nonweighted interpolation on the subintervals, which in view of (2.5) is also an admissible mesh for weighted interpolation with the same bound, the weight function being piecewise constant. Observe that the presence of a dominant weight forces the nonweighted error to be much smaller on the corresponding band than on the other bands.

EXAMPLE 3.8. (complex digitalfilters) Algorithm AFP can be easily adapted to produce polynomial filters in the complex plane. We consider the following example, taken form [26], of a nearly linear-phase low-pass filter on the union K = [K.sub.1] [union] [K.sub.2] of two disjoint arcs of the complex unit circle,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Again, we have chosen as starting mesh the union of two admissible meshes of [K.sub.1] and [K.sub.2], which in view of (2.5) is also an admissible mesh for weighted interpolation with the same bound. As in EXAMPLE 3.4, we use the standard complex monomial basis to construct the Vandermonde matrix, with two refinement iterations in Algorithm AFP.

In Figure 3 18 we show the approximate Fekete points of degree n = 31 and the estimated Lebesgue constants up to degree n = 35. The interpolation error in the weighted norm at degree [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], to be compared with an error of about 0.04 obtained in [26] with the optimal polynomial filter of degree 31, computed by the Remez algorithm.

EXAMPLE 3.9. (two-dimensional digitalfilters) Two-dimensional digital filters have important applications in the processing of images and other two-dimensional signals; see, e.g., [1, 28]. Here we show two examples of two-dimensional filters constructed by interpolation at approximate Fekete points. We begin with the following response function on a square domain,

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

with a square passband and a square frame as stopband; see Figures 3.19-3.20. These kinds of "rectangularly symmetric" two-dimensional filters are discussed, e.g., in [28].

[FIGURE 3.14 OMITTED]

[FIGURE 3.15 OMITTED]

In order to compute approximate Fekete points of not small degree by algorithm AFP, already in two dimensions it begins to be important to start from a weakly admissible mesh instead of an admissible mesh; cf. [9, 15]. Indeed, two-dimensional admissible meshes of degree n obtained by Markov inequalities have O(n4) points, due to a O(1/n2) spacing, and this leads to a heavy computational load in algorithm AFP.

[FIGURE 3.16 OMITTED]

On the other hand, weakly admissible meshes are discrete subsets of a compact K, where a polynomial inequality like (2.3) holds where C = [C.sub.n] is not constant but increases at most algebraically with n. Such a relaxation of the polynomial inequality implies that their cardinality can be much lower than that of admissible meshes. Recall that, for example, any set X of cardinality N = dim([P.sup.d.sub.n](K)), which is unisolvent for interpolation of degree n, satisfies (2.3) with C = [C.sub.n] equal to the Lebesgue constant (thus, e.g., Fekete points of K form a weakly admissible mesh). As with admissible meshes, weakly admissible meshes can be constructed by finite unions. All the relevant inequalities still hold with [C.sub.n] replacing C.

In the present example we have used the weakly admissible mesh obtained by union of the Padua points of the internal square and of four rectangles giving the square frame. We recall that the Padua points are the first known example of optimal points for total degree polynomial interpolation in two variables, with a Lebesgue constant increasing like log squared of the degree; cf. [7, 8, 14]. This implies that in (2.3>(2.5) we have C = [C.sub.n] = O([log.sup.2] n).

In Figure 3.19 we show the two-dimensional nonweighted polynomial filters obtained by interpolation at approximate Fekete points, extracted from such weakly admissible meshes at degree n = 20 and n = 30 (the interpolation errors on K are about 0.13 and 0.04, respectively). In Figure 3.20 we see the approximate Fekete points of degree n = 20 with the underlying weakly admissible mesh made of Padua points and the growth of the estimated Lebesgue constants.

[FIGURE 3.17 OMITTED]

[FIGURE 3.18 OMITTED]

In this simple example, a filter could have been constructed also as tensor-product of one-dimensional filters (cf. [1]), but the difference is that the present is a total-degree filter, which for the same n has roughly half the number of coefficients, namely N = dim([P.sup.2.sub.n](K)) = (n + 1)(n + 2)/2 instead of [(n + 1).sup.2] = dim([P.sup.1.sub.n](K) [cross product] [P.sup.1.sub.n](K)) coefficients.

[FIGURE 3.19 OMITTED]

[FIGURE 3.20 OMITTED]

A more difficult example is given by the following weighted response function,

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where the stopband and the passband are triangular. Again, we extract the approximate Fekete points from the union of weakly admissible meshes of the two triangles. In [9], the concept of "geometric" weakly admissible mesh has been developed, namely one obtained by a suitable geometric transformation from a known weakly admissible mesh on a reference domain. In the case of a triangle, we can simply map the Padua points of degree 2n in [[- 1, 1].sup.2] to the triangle by the well-known Duffy quadratic transformation [19], obtaining a weakly admissible mesh of degree n for the triangle with [C.sub.n] = O ([log.sup.2] 2n); see [9] for a more complete discussion.

In Figure 3.21 we show the corresponding nonweighted (left) and weighted ([w.sub.1] = 1, [w.sub.2] = 10; right) filters of degree n = 20. The interpolation error of the nonweighted filter is about 0.11 on the passband and 0.07 on the stopband, whereas the nonweighted error of the weighted filter becomes about 0.28 on the passband and is reduced to about 0.01 on the stopband. The respective approximate Fekete points are plotted in Figure 3.22, and the estimated Lebesgue constants in Figure 3.23.

[FIGURE 3.21 OMITTED]

[FIGURE 3.22 OMITTED]

[FIGURE 3.23 OMITTED]

* Received June 8, 2009. Accepted for publication October 22, 2009. Published online on February 28, 2010. Recommended by F. Marcellan. Supported by the "ex-60%" funds and by the biennial Project "Interpolation and Extrapolation: Theory, Algorithms and Applications" of the University of Padova, and by the INdAM-GNCS.

REFERENCES

[1] A. ANTONIOU AND W.-S. LU, Two-Dimensional Digital Filters, Marcel Dekker, New York, 1992.

[2] T. BAGBY, On interpolation by rational functions, Duke Math. J., 36 (1969), pp. 95-104.

[3] J. BAGLAMA, D. CALVETTI, AND L. REICHEL, Fast Leja points, Electron. Trans. Numer. Anal., 7 (1998), pp. 124-140. http://etna.math.kent.edu/vol.7.1998/pp124-140.dir/pp124-140.html.

[4] J. P. BERRUT, The barycentric weights of rational interpolation with prescribed poles, J. Comput. Appl. Math., 86 (1997), pp. 45-52.

[5] J. P. BERRUT, R. BALTENSPERGER, AND H. D. MITTELMANN, Recent developments in barycentric rational interpolation, Trends and Applications in Constructive Approximation, pp. 27-51, Internat. Ser. Numer. Math., 151, Birkhauser, Basel, 2005.

[6] P. BORWEIN AND T. ERDELYI, Polynomials and Polynomial Inequalities, Springer, New York, 1995.

[7] L. BOS, M. CALIARI, S. DE MARCHI, M. VIANELLO, AND Y. XU, Bivariate Lagrange interpolation at the Padua points: the generating curve approach, J. Approx. Theory, 143 (2006), pp. 15-25.

[8] L. BOS, S. DE MARCHI, M. VIANELLO, AND Y. XU, Bivariate Lagrange interpolation at the Padua points: the ideal theory approach, Numer. Math., 108 (2007), pp. 43-57.

[9] L. BOS, J.P. CALVI, N. LEVENBERG, A. SOMMARIVA, AND M. VIANELLO, Geometric weakly admissible meshes, discrete least squares approximation and approximate Fekete points, preprint, 2009. Available online at: http://www.math.unipd.it/~marcov/publications.html.

[10] L. BOS AND N. LEVENBERG, On the approximate calculation of Fekete points: the univariate case, Electron. Trans. Numer. Anal., 30 (2008), pp. 377-397. http://etna.math.kent.edu/vol.30.2008/pp377-397.dir/pp377-397.html.

[11] L. BOS, N. LEVENBERG, AND S. WALDRON, On the spacing of Fekete points for a sphere, ball or simplex, Indag. Math., 19 (2008), pp. 163-176.

[12] P. A. BUSINGER AND G. H. GOLUB, Linear least-squares solutions by Householder transformations, Numer. Math., 7 (1965), pp. 269-276.

[13] C. S. BURRUS, Digital Signal Processing and Digital Filter Design, (Draft), Sept. 14, 2009. Available at Connexions Web site: http://cnx.org/content/col10598/1.5/.

[14] M. CALIARI, S. DE MARCHI, AND M. VIANELLO, Bivariate polynomial interpolation in the square at new nodal sets, Appl. Math. Comput., 165 (2005), pp. 261-274.

[15] J. P. CALVI AND N. LEVENBERG, Uniform approximation by discrete least squares polynomials, J. Approx. Theory, 152 (2008), pp. 82-100.

[16] C. CARSTENSEN AND G. MUHLBACH, The Neville-Aitken formula for rational interpolants with prescribed poles, Numer. Algorithms, 3 (1992), pp. 133-141.

[17] A. CIVRIL AND M. MAGDON-ISMAIL, On selecting a maximum volume sub-matrix of a matrix and related problems, Theoret. Comput. Sci., 410 (2009), pp. 4801-4811.

[18] K. DECKERS, J. VAN DEUN, AND A. BULTHEEL, Rational Gauss-Chebyshev quadrature formulas for complex poles outside [-1, 1], Math. Comp., 77 (2008), pp. 967-983.

[19] M. G. DUFFY, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAMJ. Numer. Anal., 19 (1982), pp. 1260-1262.

[20] C. F. DUNKL AND Y. XU, Orthogonal Polynomials of Several Variables, Encyclopedia of Mathematics and its Applications 81, Cambridge University Press, Cambridge, 2001.

[21] M. EMBREE AND L. N. TREFETHEN, Green's functions for multiply connected domains via conformal mapping, SIAM Rev., 41 (1999), pp. 745-761.

[22] W. GAUTSCHI, The use of rational functions in numerical quadrature, J. Comput. Appl. Math., 133 (2001), pp. 111-126.

[23] M. HASSON, The degree of approximation by polynomials on some disjoint intervals in the complex plane, J. Approx. Theory, 144 (2007), pp. 119-132.

[24] M. HASSON AND J. M. RESTREPO, Approximation on disjoint intervals and its applicability to matrix preconditioning, Complex Var. Elliptic Equ., 52 (2007), pp. 757-769.

[25] LAPACK,Users' Guide, SIAM, Philadelphia, 1999.

[26] B. LE BAILLY AND J. P. THIRAN, Computing complex polynomial Chebyshev approximants on the unit circle by the real Remez algorithm, SIAM J. Numer. Anal., 36 (1999), pp. 1858-1877.

[27] A. L. LEVIN AND E. B. SAFF, Potential theoretic tools in polynomial and rational approximation, in Harmonic Analysis and Rational Approximation, J.-D. Fournier, J. Grimm, J. Leblond, and J. R. Partington, eds., Lecture Notes in Control and Inform. Sci., 327, Springer, 2006, pp. 71-94.

[28] W.-S. LU AND A. ANTONIOU, Design of Two-Dimensional Digital Filters, in CRC Handbook of Electrical Filters, J. T. Taylor and Q. Huang, eds., CRC Press, Boca Raton, FL, 1997, pp. 185-206.

[29] Y. MADAY, N. C. NGUYEN, A. T. PATERA, AND G. S. H. PAU, A general multipurpose interpolation procedure: the magic points, Commun. Pure Appl. Anal., 8 (2009), pp. 383-404.

[30] G. MASTROIANNI AND G. V. MILOVANOVIC,Weighted interpolation of functions with isolated singularities, in Approximation Theory: A Volume Dedicated to Blagovest Sendov, B. Bojanov, ed., DARBA, Sofia, 2002, pp. 310-341.

[31] THE MATHWORKS,MATLAB Documentation Set, Natick, MA, 2009 version.

[32] T. W. PARKS AND J. H. MCCLELLAN, Chebyshev approximation for nonrecursive digital filters with linear phase, IEEE Trans. Circuit Theory, Vol. CT-19, no. 2 (1972), pp. 189-194.

[33] W. PLESNIAK, Remarks on Jackson's theorem in RN, East J. Approx., 2 (1996), pp. 301-308.

[34] L. REICHEL, Some computational aspects of a method for rational approximation, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 1041-1057.

[35] Y. SAAD, Filtered conjugate residual-type algorithms with applications, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 845-870.

[36] E. B. SAFF AND V. TOTIK, Logarithmic Potentials with External Fields, Springer, Berlin, 2007.

[37] B. V. SHABAT, Introduction to Complex Analysis. Part II. Functions of Several Variables, Translations of Mathematical Monographs 110, American Mathematical Society, Providence, 1992.

[38] J. SHEN, G. STRANG, AND A. J. WATHEN, The potential theory of several intervals and its applications, Appl. Math. Optim., 44 (2001), pp. 67-85.

[39] A. SOMMARIVA AND M. VIANELLO, Computing approximate Fekete points by QR factorizations of Vandermonde matrices, Comput. Math. Appl. 57 (2009), pp. 1324-1336.

[40] A. SOMMARIVA AND M. VIANELLO, Gauss-Green cubature and moment computation over arbitrary geometries, J. Comput. Appl. Math. 231 (2009), pp. 886-896.

[41] J. VAN DEUN, Eigenvalue problems to compute almost optimal points for rational interpolation with prescribed poles, Numer. Algorithms, 45 (2007), pp. 89-99.

[42] J. L. WALSH, Interpolation and approximation by rational functions in the complex domain. Fourth ed. American Mathematical Society Colloquium Publications, Vol. XX, American Mathematical Society, Providence, 1965.

[43] J. L. WALSH AND H. G. RUSSELL, On the convergence and overconvergence of sequences of polynomials of best simultaneous approximation to several functions analytic in distinct regions, Trans. Amer. Math. Soc., 36 (1934), pp. 13-28.

[44] T. WARBURTON, An explicit construction of interpolation nodes on the simplex, J. Engrg. Math., 56 (2006), pp. 247-262.

[45] P. ZITNAN, A stable collocation solution of the Poisson problems on planar domains, Abstracts of the 2nd Dolomites Workshop on Constructive Approximation and Applications, Canazei (Italy), Sept. 2009. Available online at: http://www.math.unipd.it/~dwcaa09

A. SOMMARIVA ([dagger]) and M. VIANELLO ([dagger])

([dagger]) Department of Pure and Applied Mathematics University of Padova, via Trieste 63, 35121 Padova, Italy ({alvise, marcov}@math.unipd.it).
COPYRIGHT 2010 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.