On the calculation of approximate Fekete points: the univariate case.

AMS subject classifications. 41A05, 42A15, 65D05

1. Introduction. Fekete points are a set of points which are good for polynomial interpolation that may be defined for any compact set in any dimension. They are therefore a natural and important set of points from the point of view of polynomial interpolation. How ever, their computation involves an expensive multivariate optimization which is moreover numerically challenging for higher degrees; cf. the webpage of Womersley . Recently Sommariva and Vianello  proposed a method for calculating approximate Fekete points that is highly efficient and may easily be used on different compact sets.

Polynomial interpolation, at least in one variable, is a classical subject. However, to make the notions we consider here precise, we briefly outline the main features of (multivariate) polynomial interpolation. Consider K [subset] [R.sup.d] a compact set. The polynomials of degree at most n in d real variables, when restricted to K, form a certain vector space which we will denote by [P.sub.n] (K). The space [P.sub.n] (K) has a dimension [N.sub.n] := dim([P.sub.n] (K)). The polynomial interpolation problem for K is then, given a set of [N.sub.n] distinct points [A.sub.n] [subset] K and a function f : K [right arrow] R, to find a polynomial p [member of] [P.sub.n] (K) such that

(1.1) p(a) = f(a), [for all]a [member of] [A.sub.n].

If we choose a basis,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

of [P.sub.n] (K), then any polynomial p [member of] [P.sub.n] (K) may be written in the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for some constants [c.sub.j] [member of] R. Hence the conditions (1.1) may be expressed as

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which are exactly [N.sub.n] linear equations in [N.sub.n] unknowns [c.sub.j]. In matrix form this becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where c [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the vector formed of the [c.sub.j] and F is the vector of function values f(a), a [member of] [A.sub.n] This linear system has a unique solution precisely when the so-called Vandermonde determinant

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If this is the case, then the interpolation problem (1.1) is said to be correct.

Supposing then that the interpolation problem (1.1) is correct, we may write the interpolating polynomial in so-called Lagrange form as follows. For a [member of] [A.sub.n] set

(1.4) [l.sub.a](x) := vdm([A.sud.n]\{a}[union]{x};[B.sub.n])/vdm([A.sub.n];[B.sub.n])

Note that the numerator is just the Vandermonde determinant with the interpolation point a [member of] [A.sub.n] replaced by the variable x [member of] [R.sup.d].

It is easy to see that [l.sub.a] [member of] [P.sub.n] (K). Moreover [l.sub.a](b) = [[delta].sub.ab], the Kronecker delta, for b [member of] [A.sub.n]. Using these so-called Fundamental Lagrange Interpolating Polynomials we may write the interpolant of (1.1) as

(1.5) p(x) [summation over a[member of][A.sub.n] f(a)[l.sub.a](x).

The mapping f [right arrow] p is a projection and hence we write [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. If we regard both f, p [member of] C(K) then the operator [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] has operator norm (as is not difficult to see)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This operator norm, called the Lebesgue constant, gives a bound on how far the interpolant is from the best uniform polynomial approximant to f. It follows that the quality of approximation to f provided by the interpolant p is indicated by the size of the Lebsegue constant, the smaller it is the better.

Now, suppose that [F.sub.n] [subset] K is a subset of [N.sub.n] distinct points for which [A.sub.n] = [F.sub.n] maximizes [absolute value of vdm([A.sub.n]; [B.sub.n])]. Then by (1.4), each supremum norm

(1.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and hence the corresponding Lebesgue constants are such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

i.e., the Lebesgue constants grow polynomially in n, which is the best that is known in general. Such a set [F.sub.n] (it may not be unique) is called a set of (true) Fekete points of degree n for K and provide, for any K, a good (typically excellent) set of interpolation points. For more on Fekete points (and polynomial interpolation) we refer the reader to  (and its references). Note that the Fekete point sets [F.sub.n] and also the Lebesgue constants [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are independent of the basis [B.sub.n]. We also remark that for each degree n, the Fekete points [F.sub.n], form a set, i.e., they do not provide an ordering of the points. In applications, especially for high degrees, the ordering of the points can be important. One such ordering is provided by the so-called Leja points; see, for example, [10, 1]. Note, however, that the Leja points provide an ordering of all the points up to those of degree n, whereas the Fekete points are for only degree n. Once a set of Fekete points has been calculated they can be ordered by the Leja method, but we do not pursue that here.

Now for the Sommariva-Vianello method. Here is some MATLAB code that implements the algorithm to find 21 (degree 20) approximate Fekete interpolation points for the interval [-1,1]:

n=21; % number of interpolation points m=1000; x=linspace(-1,1,m); % discrete model of [-1,1] A = gallery('chebvand',n,x); % A is the n by m Vandermonde matrix % in the Chebyshev polynomial basis b=rand(n,1); % a random rhs y=A\b; % y is the MATLAB solution of Ay=b pp=y~=0; % vector of indices of the non-zero elements of y pts=x(pp) % selects the points from x according to pp

Before explaining what this code does, we give some example results. In Figure 1.1 we show the approximate Fekete points computed by the above code (marked by +) as well as the true Fekete points, the extended Chebyshev points and also the so-called Leja points. Note how remarkably close the approximate Fekete points are to the true Fekete points. In comparison, the Vandermonde determinant (using the Chebyshev basis) for the true Fekete points is 1.532. [10.sup.11] (the maximum possible), whereas, for the approximate Fekete points, it is 1.503. [10.sup.11]. For the extended Chebyshev points (a well-known set of excellent interpolation points for [-1,1]), this determinant is 1.265. [10.sup.11]. As we shall see in Theorem 4.1 below, the values of the Vandermonde determinants for the approximate Fekete points are sufficiently close to the values of the Vandermonde determinants for the true Fekete points to allow us to conclude that these two sets of points are asymptotically the same.

Further, the Lebesgue constant for the true Fekete points is approximately 2.6, while for the approximate Fekete points it is approximately 2.8 and for the extended Chebyshev points approximately 2.9. The reader may be interested to note that for 21 equally spaced points on the interval [-1,1] the Lebesgue constant is approximately 10986.5, i.e., significantly greater. We note the classical results that for n+1 Chebyshev points in the interval and 2n+1 equally spaced points on a circle, the Lebesgue constants grow like c In (n) (and this order of growth is optimal). This is further discussed in Section 3.2.

We hope that these results convince the reader that the Sommariva-Vianello algorithm is indeed promising. The purpose of this paper is to discuss some of the theoretical aspects of the algorithm; to hopefully explain why this algorithm gives such good results. Numerical implementation is discussed in . In the next section we show how the procedure is related to a natural greedy algorithm for constructing submatrices of maximal determinant. We give concrete examples of the algorithm in Section 3. Finally, in Section 4 we prove that the algorithm produces points which asymptotically exhibit the same behavior as that of the true Fekete points for a finite union of nondegenerate compact, connected sets in the complex plane.

2. The relation to a greedy algorithm for maximum volume submatrices. The key to understanding how the code segment works is the command y=A\b. First observe that the command A=chebvand (n, x) produces the Vandermonde matrix of the first n Chebyshev polynomials evaluated at the points of the vector x. Specifically, the (i, j) entry of A is [T.sub.i-1]([x.sub.j]) so that the ith row of A corresponds to the Chebyshev basis polynomial [T.sub.i-1] and the jth column of A corresponds to the jth point in x, [x.sub.j]. Hence selecting a subset of columns of A corresponds to selecting a subset of the points of the vector x.

Now, note that the matrix A [member of] [R.sup.nxm] with n = 21 and m = 1000, in this case, so that the linear system Ay = b is severely underdetermined. MATLAB resolves this problem by first computing the QR factorization of A (Q [member of] [R.sup.nxn], orthogonal and R [member of] [R.sup.nxm], upper triangular) with column pivoting; cf. [7, [section]5.4]. The basics of this algorithm are easy to describe. In fact, the pivoting is a procedure to select the n "most significant" among the m > > n columns of A.

[FIGURE 1.1 OMITTED]

The first column selected is the one of largest euclidean length--call this column [a.sub.1] [member of] [R.sup.n]. Then, an orthogonal matrix [Q.sub.1] [member of] [R.sup.nxn] is chosen that maps [a.sub.l] to the first column of an upper triangular matrix, i.e., [+ or -] [parallel][a.sub.l][parallel][2.sup.e]l. Here we use the notation [parallel]x[[parallel].sub.2] to denote the euclidean ([l.sub.2]) norm of a vector in [R.sup.m] or [C.sup.m] (m may vary). Also, we use a dot "." to denote the euclidean inner product. We then compute

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [A.sub.1] [member of][R.sup.(n-1)x(m-1)]. Then these two operations are repeated to the matrix [A.sub.l], [A.sub.2] [member of] [R.sup.(n-2)x(m-2)] and so on. After n steps, we arrive at Q = [Q.sub.1] ... [Q.sub.n-2] [Q.sub.n-1] and R = [A.sub.n-1]. Once these have been calculated, MATLAB solves the system

[??][??] = b,

where [??] [member of] [R.sup.nxn] consists of the n columns so selected and [??] [member of] [R.sup.n]. The other entries of y [member of] [R.sup.m] are set to 0. Hence the command pp=y~=0; in the code segment above gives the indices of the selected columns.

Let us look now at the second column of A that the algorithm selects (the first is the one of longest euclidean length). This would correspond to the column of [A.sub.1] (in [R.sup.n-1]) of longest length. Let us call this longest column of [A.sub.1], [b.sub.1] [member of] [R.sup.n-1], and the corresponding column of A, [a.sub.2]. Then by (2.1) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for [alpha] = [+ or -] [parallel][a.sub.1][[parallel].sub.2] and some [beta] [member of] R. We have then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

from which we see that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is the component of [Q.sub.1][a.sub.2] orthogonal to [Q.sub.1][a.sub.1], and hence [parallel][b.sub.1][[parallel].sub.2] x [parallel][Q.sub.1][a.sub.1][[parallel].sub.2] is the area of the parallelogram generated by [Q.sub.1][a.sub.1], and [Q.sub.1][a.sub.2]. But [Q.sub.1] is an orthogonal matrix and so [parallel][b.sub.1][[parallel].sub.2] x [parallel][Q.sub.1][a.sub.1][[parallel].sub.2] is also the area of the parallelogram generated by [a.sub.1] and [a.sub.2]. It follows that the second column is chosen so that the area it generates with fixed al is maximal. Similarly, the third column is chosen so that the volume it generates with fixed al and [a.sub.2] is maximal, and so on.

In summary, the pivoting procedure for the QR factorization selects columns as follows:

(1) [a.sub.1] is the column of A of maximum euclidean length.

(2) Given [a.sub.1], ..., [a.sub.k], the (k + 1)st column [a.sub.k+l] is chosen so that the volume of the "box" generated by [a.sub.l], ..., [a.sub.k], [a.sub.k+l] is as large as possible.

This is precisely the standard Greedy Algorithm for constructing a maximal volume box from a collection of vectors. Note that, in principle, the algorithm selects the columns in a certain order. However, in the final result of the MATLAB command A\b this information is lost. As mentioned in the Introduction, in numerical applications the ordering of the points is important and a version of the command that saves the order information would be very useful. In this paper, however, we concentrate on the theoretical aspects of the set of points that the algorithm selects.

3. The continuous version of the algorithm. Suppose that K [subset] [R.sup.d] is compact. Given a basis [B.sub.n] = {[P.sub.1], [P.sub.2], ..., [P.sub.N]} for [P.sub.n] (K) the columns of the associated Vandermonde

matrix are of the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for some x [member of] K. Recalling that, for a Vandemonde matrix, selecting a subset of columns is equivalent to selecting a subset of points, we describe a continuous version of the Sommariva-Vianello Algorithm as follows:

(1) The first point [x.sub.i] [member of] K is chosen to maximize [parallel][??](x)[[parallel].sub.2].

(2) Given [x.sub.i], [x.sub.2], ..., [x.sub.k] the (k + 1)st point [x.sub.k+l] [member of] K is chosen so that the volume generated by the columns [??]([x.sub.k+1]) and [??]([x.sub.1]), [??]([x.sub.2]), ..., [??]([x.sub.k]) is as large as possible.

3.1. First example: the unit circle. Here we take K to be the unit circle so that [P.sub.n] (K) is the trigonometric polynomials of degree n with dimension N = 2n + 1. We take the orthonormal basis,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with respect to the inner-product,

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In particular, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

More generally,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which the reader will notice is the reproducing kernel for the space of trigonometric polynomials of degree at most n, equipped with the inner product (3.1).

It follows that for any set of 2n + 1 equally spaced angles {[[theta].sub.k]}, k = 0, 1, ..., 2n, i.e., with [[theta].sub.k] = [[theta].sub.0] + 2k[pi]/(2n + 1),

(3.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, the column vectors [??]([[theta].sub.j]) and [??]([[theta].sub.k]) are orthogonal for j [not equal to] k.

Conversely, for [theta] [member of] [0, 2[pi]] and fixed j,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

What is the result of the Algorithm in this case? Since by (3.2), the length of [??]([theta]) is constant for all [theta][member of] [0, 2[pi]], the first point chosen will be any [[theta].sub.0] [member of] [0, 2[pi]]. The second point [theta] will be so that the area generated by [??]([[theta].sub.0]) and [??]([theta]) is as large as possible. But as shown above, [??]([theta])[perpendicular to][??]([[theta].sub.0]) iff [theta] = [[theta].sub.j] for some j [not equal to] 0. Hence (noting that the lengths of [??]([theta]) are the same for all [theta][member of] R) this area will be maximized by any [[theta].sub.j], j [not equal to] 0. Continuing, we see that the output of the Algorithm is a set of equally spaced angles {[[theta].sub.j]}, generated in a random order, for some [[theta].sub.0] [member of] [0, 2[pi]]. This is also a set of true Fekete points and we see that, in this case, the approximate Fekete points of the Algorithm are even true Fekete points.

3.2. Second example: the interval [-1,1] [subset] [R.sup.1]. This example is a bit indirect. We construct a good set of interpolation points on [-1,1] by first calculating approximate Fekete points on the unit circle [S.sup.1] [subset] [R.sup.2] and projecting down, i.e., (x, y) [right arrow] x [member of] [-1,1]. Such indirect procedures are likley to be very useful also in several variables.

Since (x, [+ or -]y) project to the same point x, in order to obtain n + 1 points in [-1,1] it is natural to look for an even number, 2n, points, symmetric under the mapping y [right arrow] -y, on the unit circle.

In this case the number of points (on the unit circle), 2n, is not the dimension of trigonometric polynomials of some total degree and hence choosing the "best" basis is a bit more subtle. Here we use the basis

(3.4) [B.sub.n] := {1/[square root of 2], cos([theta]), sin([theta]), ..., cos((n - 1)[theta]), sin((n - 1)[theta]), 1/[square root of 2]cos(n[theta])}.

Some words of explanation are in order. First of all, [B.sub.n] spans the trigonometric polynomials of degree at most n with sin(n[theta]) removed. This choice is made in anticipation that the 2n equally spaced angles [[theta].sub.k] = k[pi]/n, k = 0, 1, ..., 2n-1 are the best possible for the purposes of interpolation. Note that sin(n[theta]) is equal to 0 precisely at these points, and hence must be removed from any candidate basis.

The normalization constants 1/[square root of 2] for the first and last elements of [B.sub.n] are in order that [B.sub.n] is orthonormal with respect to the equally weighted discrete inner-product based on the equally spaced points [[theta].sub.k]. Specifically, we have:

LEMMA 3.1. Suppose that [[theta].sub.k] := k[pi]/n, 0 [less than or equal to] k [less than or equal to] 2n - 1. Then the elements of [B.sub.n] are orthonormal with respect to the inner-product

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(and corresponding norm [parallel]p[parallel] = [square root of <p,p>).

Proof. We first calculate

(3.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now, for integer t, 0 < t < 2n,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Hence, by (3.6), for 0 [less than or equal to] j < m [less than or equal to] n (and consequently 0 < j + m < 2n and also 0 < [absolute value of j - m] < 2n), we have

(cos(j[theta]), cos(M[theta])) = 0.

Also, for 0 < j = m < n, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For j = m, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For j = m = 0, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The calculations with the sines are similar and so we omit them. []

With the basis [B.sub.n] the columns of the Vandermonde matrix are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that by Lemma 3.1, the rows of the matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are orthogonal and have the same euclidean length. It follows that V is a constant multiple of an orthognal matrix so that the column vectors are also orthogonal. In other words,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note also that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which is maximized precisely at [theta] = [[theta].sub.k] for some k. Moreover,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

From these considerations it follows that the Algorithm, applied to the basis [B.sub.n] and [theta] [member of] [0, 2[pi]] will select the angles [[theta].sub.k], k = 0, 1, ..., 2n - 1, in some random order.

The projected points are then cos(k[pi]/n), 0 [less than or equal to] k [less than or equal to] n, which are precisely the so-called extended Chebyshev points. Hence the Algorithm again produces an excellent set of interpolation points.

Actually, from this setup it is easy to see why the Chebyshev points have small Lebesgue constant. Although this is somewhat tangential to our main purpose it is useful to record this here as it may also shed some light on the multivariate case.

The basic fact to note is that the trigonometric (span of [B.sub.n]) fundamental Lagrange polynomial for [[theta].sub.k] is just

(3.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [K.sub.n] ([theta], [phi]) is the reproducing kernel for the span of [B.sub.n] and inner-product (3.5). That this holds is an example of a more general phenomenon. Indeed, we have for any p [member of] span([B.sub.n]),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then since p [member of] span([B.sub.n]) was arbitrary, it follows that

1/n [K.sub.n]([[theta].sub.j], [[theta].sub.k]) = [[delta].sub.jk]

and hence (1/n)[K.sub.n]([theta],[[theta].sub.k]) [member of] span([B.sub.n]) indeed coincides with the fundamental Lagrange polynomial [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] ([theta]).

We may easily compute [K.sub.n]. Indeed, by Lemma 3.1, [B.sub.n] is orthonormal with respect to the inner-product (3.5) and we have

(3.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now note that from this formula for [K.sub.n], together with (3.7) and the fact that sin(n[[theta].sub.k]) = 0, [for all]k, it follows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We may simplify

(3.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by a standard calculation.

From this it is easy to calculate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

from which it follows that the Lebesgue function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A more refined, but standard, calculation shows that actually

(3.10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for some constant c > 0, whose exact value is not important here; see, e.g., [11, [section]1.3] for similar calculations.

Further, since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We may use these symmetry properties to obtain a relation between the trigonometric fundamental Lagrange polynomials for the angles [[theta].sub.k] and the algebraic fundamental Lagrange polynomials for the points [x.sub.k] = cos([[theta].sub.k]) [member of] [-1,1]. (We emphasize that there are exactly n+1 different [x.sub.k] as [x.sub.2n-k] = [x.sub.k], as is easily seen.)

Specifically, note that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a combination of cosines only and hence even in [theta]. Therefore the substitution x = cos([theta]) results in an algebraic polynomial of degree n,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is easy to check that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and hence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the algebraic Lagrange polynomial for [x.sub.0] among the n + 1 points [x.sub.k].

Similarly, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is also even in [theta] and so

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is the algebraic fundamental Lagrange polynomial corresponding to [x.sub.n].

The construction for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is slightly different. Although in this case [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is not even, the combination

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is even. It is easy then to check that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for 0<k<n.

It follows that the algebraic Lebesgue function (for the extended Chebyshev points)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3.3. Using a general orthonormal basis. Suppose that p is a (regular) measure on [-1,1] with associated inner-product

(3.11)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We consider an orthonormal basis

[B.sub.n] := {[P.sub.0],[P.sub.1], ..., [P.sub.n] [subset] [P.sub.n]}([-1,1])

for the polynomials of degree at most n, [P.sub.n]([-1,1]).

The columns of the Vandermonde matrix are then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with norm (squared)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that this just the diagonal of the reproducing kernel for the space [P.sub.n]([-1,1]) with inner-product (3.11), sometimes also referred to as the reciprocal of the associated Christoffel function. To emphasize this relation we set

[K.sub.n](x) := [n.summation over (k=0)][P.sup.2.sub.k](X),

so that

[parallel][??](x)[parallel] = [square root of [K.sub.n](x).

To the measure [mu] is associated the corresponding Gauss-Christoffel quadrature formula. Specifically, if [a.sub.0], [a.sub.1] ..., [a.sub.n [member of] (-1,1) are the zeros of [P.sub.n+1] (x), then, for all polynomials Q(x) of degree at most 2n + 1,

(3.12) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the weights are given by

[W.sub.k] = 1/[K.sub.n]([a.sub.k])

Now consider the normalized columns

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Specializing to x = [a.sub.k], we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Define the matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The ith row of U is the vector

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that the euclidean inner-product of [U.sub.i] and [U.sub.j], by Gauss-Christoffel quadrature (3.12), is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since the polynomials Pi are orthonormal.

In other words, the rows of the matrix U are orthonormal vectors. It follows that U is an orthogonal matrix and that also the columns of U are orthonormal vectors, i.e.,

[??]([a.sub.i]) x [??]([a.sub.j]) = [[delta].sub.ij].

Now, if we apply the Sommariva-Vianello algorithm to the normalized colums, [??](x), all of length 1, the first point chosen will be random, as there is no way for the Algorithm to distinguish first points. However, if by some means, we set the first point to [x.sub.0] = [a.sub.j] for some j, 0 [less than or equal to] j [less than or equal to] n, then since the vectors [??]([a.sub.i]), i [not equal to] j, are orthogonal to [??]([a.sub.j]) (and to each other), the remaining n points selected by the algorithm will be just the [a.sub.i], i [not equal to] j, in some random order.

To summarize, if the Sommariva-Vianello algorithm, using an orthonormal basis, is applied to normalized columns and the first point is (artificially) set to one of the Gauss-Christoffel quadrature points, then the Algorithm selects precisely the Gauss-Christoffel quadrature points associated to the measure of orthogonality. We remark that if, on the other hand, the first point is not so set, the first point would be randomly chosen in the interval and would in general be a poor choice (although the algorithm would eventually self-correct).

3.3.1. An illustrative example. When we normalize the columns to have length one, it is also possible that the Algorithm selects good interpolation points, other than the Gauss-Christoffel points, depending on how the first point is set. For example, consider the measure

d[mu] = 1/[square root of 1 - [x.sup.2]] dx.

The orthonormal polynomials with respect to this measure are the Chebyshev polynomials of the first kind. Specifically,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If we set x = cos([theta]) and y = cos([phi]), then we may compute

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by a standard calculation.

We note that the Gauss-Christoffel quadrature points are the zeros of [T.sub.n+1](x), i.e.,

[a.sub.k] = cos(2k + 1/2(n + 1)[pi]), k = 0, 1, ..., n.

By the general theory of the previous section (or else by direct calculation), it follows that {[??]([a.sub.k])} is a set of orthogonal vectors.

But notice that by (3.13), if we set

[b.sub.k] = cos(2k/2n + 1 [pi]), k = 0,1, ..., n,

then the set of vectors

{[??]([[theta].sub.k]) : k = 0,1, ...,n}

is also orthogonal. It follows that if we (artificially) set the first point to [x.sub.0] = [b.sub.0] = cos(0) = 1, a not unnatural choice, then the Algorithm applied to normalized colums will select the points [b.sub.k] in some random order.

Note, however, that these selected points are not symmetric, specificially, x = 1 is included, but x = -1 is not. We record, for comparison's sake, that for 21 points (degree 20), the Vandermonde determinant is approximately 5.796 - [10.sup.10] and the Lebesgue constant approximately 3.3.

3.4. Gram determinants and other bases. The volume of the parallelotope generated by the vectors [[??].sub.1], [[??].sub.2], ..., [[??].sub.k] [member of] RI is given by the associated Gram determinant. Specifically,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For more information on Gram determinants, see, for example, [6, [section]8.7].

The continuous version of the Sommariva-Vianello algorithm can then be re-phrased as

(1) The first point [x.sub.1] [member of] K is chosen to maximize [parallel][??](x)][[parallel].sub.2].

(2) Given [x.sub.1], [x.sub.2], ..., [x.sub.k] the (k + 1)st point [x.sub.k+1] [member of] K is chosen so that the Gram determinant,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is as large as possible.

3.4.1. The standard basis. In this section we briefly consider the use of the standard polynomial basis,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In this case the columns of the Vandermonde matrix are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then we compute

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In particular, for x = y,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Clearly, [parallel][??](x)[parallel] is maximized for x = [+ or -]1 and the algorithm will choose one of these at random. Let us suppose without loss of generality that the first point chosen is [x.sub.1] = +1. In general, it is difficult to calculate the precise points that the algorithm chooses. For n even it is not too difficult to show that the second point is [x.sub.2] = -1. Specifically, the second point will be the one for which g([??](1),[??](x)) is maximal. But we calculate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We claim that, for n even, this is maximized for x = -1. To see this, first note that then

[??](-1) x [??](1) = [(-1).sup.n] - 1/-1-1 = 0

and hence, for all x [member of] [-1,1],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, g([??](1),[??](x)) is indeed maximized for x = -1 and [x.sub.2] = -1 will be the second point chosen.

If n is odd, then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and already proving that [x.sub.2] = -1 (although we conjecture that is the case) seems to be a bit difficult.

To see the effect of the basis, we will compute the four points chosen for n = 4 using both the standard and Chebyshev basis.

Continuing with the standard basis, since n = 4 is even we already know that [x.sub.1] = +1 (by choice) and [x.sub.2] = -1. Let us compute [x.sub.3]. This will be the point for which g([??](1),[??](-1), [??](x)) is maximized. But

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The derivative of this expression is

d/dx 8(1 - [x.sup.4])(1 - [x.sup.2]) = 16x ([X.sup.2] - 1)(3[x.sup.2] + 1),

and hence the Gram determinant is maximized at x = 0 and [x.sub.3] = 0.

To find the fourth point, we must maximize

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The derivative of this determinant is

d/dx 4[x.sup.2][([x.sup.2] - 1).sup.2] = 8x([x.sup.2] - 1)(3[x.sup.2] - 1).

Hence the fourth point is [x.sub.4] = [+ or -]1/[square root of 3]. In summary, the four points for n = 4 and the standard basis are

{1,-1,0,[+ or -] [+ or -]/[square root of 3]}

which are not particularly good for interpolation.

Now let us compute the points for n = 4 using the Chebyshev basis

[B.sub.4] = {[T.sub.0](x), [T.sub.1](x), [T.sub.2](x), [T.sub.3](x)},

so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is easy to check that again [x.sub.1] = 1 (by choice) and [x.sub.2] = -1. To find the third point [x.sub.3], we must maximize

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The derivative of this determinant is

d/dx 32[([x.sup.2] - 1).sup.2](4[x.sup.2] + 1) = 128x ([x.sup.2] - 1)(6[x.sup.2] - 1).

Hence the third point will be [x.sub.3] = [+ or -]1/[square root of]- = .4082482...

If we take [x.sub.3] = +1/[square root of 6], then, after some tedious calculations, we find that the fourth point [x.sub.4] = ([square root of 6] - [square root of 114])/18 = -.457088...

For comparison's sake, the Fekete points are

{-1, 1/[square root of 5], + 1/[square root of 5] + {-1, -.4472135954, +.4472135954, +1},

and we see that using the Chebyshev basis has produced a better result than for the standard basis.

We mention in passing that for the interval there is always an artificial basis for which the Algorithm selects precisely the true Fekete points. This is, of course, not very useful as one must know the true Fekete points a priori to construct the basis.

PROPOSITION 3.2. Suppose that K = [-1,1] and that [B.sub.n] is the Lagrange basis based at the true Fekete points for K. Then the Sommariva-Vianello algorithm will select the true Fekete points in some order.

Proof. Let us denote the set of true Fekete points by

[F.sub.n] = {[a.sub.0], [a.sub.1], ..., [a.sub.n]}.

Then,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that [??]([a.sub.j]) = [[??].sub.j], the standard basis vector. Hence the vectors [??]([a.sub.j]) are mutually orthogonal. Moreover, by Fejer's Theorem (cf. ),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [parallel][??](x)[[parallel].sub.2] = 1 only for x [member of] [F.sub.n] The result folows. []

We remark that the same result holds in several variables for any compact set K and degree n with the property that the sum of the squares of the fundamental Lagrange polynomials based at the Fekete points is at most one, and attains 1 only at the Fekete points (an unfortunately rare occurence; see, e.g.,  for a discussion of this problem).

Of course, in applications, it will be important to select the "right" basis. In general, we would suggest the use of a basis of polynomials orthonormal with respect to the so-called equilibrium measure for K of potential theory. (For K = [-1,1] this would be 1/[pi] 1/[square root of 1 - [x.sup.2]] dx and the orthonormal polynomials are (multiples) of the Chebyshev polynomials; for K the unit circle this measure is just 1/2[pi]d[theta] and the orthonormal polynomials are (multiples of) the trigonometric polynomials.) However, for general K this measure (and the associated polynomials) are difficult to calculate. Sommariva and Vianello  discuss a form of "iterative improvement" that adjusts the basis dynamically.

4. A convergence theorem for K [subset] C. In this section, we work in the complex plane (C. Although the basis used does have a strong influence on the points the Algorithm selects, it turns out that asymptotically the points will have the same distribution, that of the true Fekete points.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

denote the Vandermonde determinant for these points with respect to the standard basis. For K [subset] C compact, we call

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the transfinite diameter of K. Thus if {[f.sub.1.sup.(n)], [f.sub.2.sup.(n)], ..., [f.sub.n+1.sup.(n)]} is a set of true Fekete points of degree n for K,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We formulate a general convergence theorem.

THEOREM 4.1. Suppose that K = [U.sub.i][K.sub.i] [subset] C is a finite union of continua (i.e., each [K.sub.i] is compact and connected, not a single point). Suppose further that [phi] : [Z.sub.+] [right arrow] R has the property that

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and that [A.sub.n] [subset] K, n = 1, 2, ..., are discrete subsets of K such that for all x [member of] K,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then for any basis the points [b.sub.l], ..., [b.sub.n+l] [member of] [A.sub.n] [subset] K generated by the Sommariva-Vianello Algorithm based on points in An satisfy

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the discrete probability measures [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] converge [weak-.sup.*] to the potential-theoretic equilibrium measure d[mu]K of K.

REMARK 4.2. For K = [-1,1], d[[mu].sub.[-1,1]] = 1/[square root of 1 - [x.sup.z]] dx; for K the unit circle [S.sup.1], d[mu][S.sup.1] = 1/2[pi]d[theta]. We refer the reader to  for more on complex potential theory. As an example of the condition (4.1), if K = [-1,1], An consisting of order [n.sup.2+[epsilon]], [member of] > 0, equally spaced points would suffice.

Proof. Suppose that {[f.sub.1.sup.(n)], [f.sub.2.sup.(n)], ..., [f.sub.n+1.sup.(n)]} is a set of true Fekete points of degree n for K. Let then {[a.sub.l], [a.sub.2], ..., [a.sub.n+1]} [subset] [A.sub.n] be such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the existence of which is guaranteed by our hypotheses.

By a result of Kovari and Pommerenke  it follows that there is a constant [c.sub.1] > 0 such that

(4.2) [absolute value of [f.sub.i.sup.(n)] - [f.sub.j.sup.(n)]] [greater than or equal to] [C.sub.1]/[n.sup.2], i [not equal to] j.

Consequently, for i [not equal to] j,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since [phi](n) = [??](1/[n.sup.2]) by (4.1).

We will first show that

(4.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

To see this, first note that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by the definition of Fekete points.

Also,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we have set [c.sub.3] := 2/[c.sub.2] > 0.

Hence, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and then (4.4) follows by taking [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] roots, and using (4.1).

Continuing, suppose that the Sommariva-Vianello Algorithm selects the points [b.sub.1], [b.sub.2] ..., [b.sub.n+l] [member of] [A.sub.n] [subset] K. Then Civril and Magdon-Ismail  have shown that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [a.sup.*.sub.1], ..., [a.sup.*.sub.n+1] are the points among [A.sub.n] which maximize the Vandermonde determinant. Note that such an inequality is independent of the basis used in calculating the determinant. Hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by the definition of the [a.sup.*.sub.j]. Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The final statement, that [[mu].sub.n] converges [weak-.sup.*] to d[[mu].sub.K], then follows by [2, Theorem 1.5].

* Received May 2, 2008. Accepted for publication October 20, 2008. Published online on December 17, 2008. Recommended by L. Reichel.

REFERENCES

 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.

 T. BLOOM, L. BOS, C. CHRISTENSEN, AND N. LEVENBERG, Polynomial interpolation of holomorphic functions in C and [C.sup.n], Rocky Mountain J. Math., 22 (1992), pp. 441-470.

 L. BOS, Some remarks on the Fejer problem for Lagrange interpolation in several variables, J. Approx. Theory, 60 (1990), pp. 133-140.

 L. BOS, N. LEVENBERG, AND S. WALDRON, On the spacing of Fekete points for a sphere, ball, or simplex, to appear in Indag. Math. (N. S.).

 A. CIVRIL AND M. MAGDON-ISMAIL, Finding maximum Volume sub-matrices of a matrix, Technical Report 07-08, Department of Computer Science, Rensselaer Polytechnic Institute, Ithaca, NY, 2007.

 P. J. DAVIS, Interpolation and Approximation, Dover, New York, 1975.

 G. GOLUB AND C. F. VAN LOAN, Matrix Computations, Third ed., Johns Hopkins U. Press, Baltimore, MD, 1996.

 T. KOVARI AND C. POMMERENKE On the distribution of Fekete points, Mathematika, 15 (1968), pp. 70-75.

 T. RANSFORD, Potential Theory in the Complex Plane, Cambridge University Press, Cambridge, 1995.

 L. REICHEL, Newton interpolation at Leja points, BIT, 30 (1990), pp. 332-346.

 T. RIVLIN, Chebyshev Polynomials, Second ed., Wiley, New York, 1990.

 A. SOMMARIVA AND M. VIANELLO, Computing approximate Fekete points by QR factoriziations of Vandermonde matrices, Comput. Math. Appl., to appear.

 R. WOMERSLEY, http://web.maths.unsw.edu.au/_rsw/Sphere/.

L. P. BOS ([dagger]) AND N. LEVENBERG ([double dagger])

([dagger]) Department of Mathematics and Statistics, University of Calgary, Calgary, Alberta, Canada T2N 1N4 (1pbos@uca1gary.ca).

([double dagger]) Department of Mathematics, Indiana University, Bloomington, Indiana, USA (nlevenbe@indiana.edu).
COPYRIGHT 2008 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.