Modified Adaptive Algorithms.
It is well known that the adaptive algorithm is simple and easy to program but the results are not fully competitive with other nonlinear methods such as free knot spline approximation. We modify the algorithm to take full advantages of nonlinear approximation. The new algorithms have the same approximation order as other nonlinear methods, which is proved by characterizing their approximation spaces. One of our algorithms is implemented on the computer, with numerical results illustrated by figures and tables.
Introduction
. It is common knowledge that nonlinear approximation methods
are better, in general, than their linear counterparts. In the case of splines,
nonlinear approximation puts more knots where the function to be approximated
changes rapidly, which results in dramatic improvements in approximating functions
with singularities. There are various satisfactory results on free knot spline approxi-
mation, in which knots are chosen at one's will. Most related theorems are proved by
showing the existence of certain balanced partitions (a more accurate description will
be given later). This may cause di#culties in practice, since it is often numerically
expensive to find such balanced partitions. Then, there is so-called adaptive approximation
by piecewise polynomial (PP) functions, in which only dyadic intervals are
used in the partition. Adaptive approximation draws great attention because of its
simplicity in nature. As a price to pay for the simplicity, its approximation power
is slightly lower than that of its free knot counterpart. Moreover, it is not known
exactly what kind of functions can be approximated to a prescribed order; that is,
there is no characterization of adaptive approximation spaces. We point out here that
when we say adaptive algorithms in this paper, we mean those that approximate a
given (univariate) function by PP functions/splines. There are other kinds of adaptive
algorithms; some are characterized in the literature (see [10] for an example).
In this paper, we shall modify the existing adaptive algorithms in two ways.
The resulting algorithms have the same approximation power as free knot spline
approximation while largely keeping the simplicity of adaptive approximation. In the
next section, we shall state some known results on free knot spline approximation.
After describing our algorithms in section 3, in section 4 we shall give our main results,
which are parallel to those on free knot spline approximation given in the next section.
Numerical implementation and examples will be the contents of the last section.
# Received by the editors March 17, 1999; accepted for publication (in revised form) February 9,
2000; published electronically August 29, 2000.
http://www.siam.org/journals/sinum/38-3/35356.html
Department of Mathematics and Computer Science, Georgia Southern University, Statesboro,
GA 30460 (yhu@gasou.edu).
# Department of Mathematics, University of Manitoba, Winnipeg, MB, Canada R3T 2N2
(kkopotun@math.vanderbilt.edu). This author was supported by NSF grant DMS 9705638.
- Department of Mathematics, Southwest Missouri State University, Springfield, MO 65804
(xmy944f@mail.smsu.edu).
1014 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
We emphasize that we consider only the univariate case in this paper. The idea
of merging cubes was initially introduced and used by Cohen et al. in their recent
paper [11] on multivariate adaptive approximation. The resulting partition consists
of rings, which are cubes with a (possibly empty) subcube removed. Their algorithm
produces near minimizers in extremal problems related to the space BV (R 2 ). The
authors further explored this algorithm in [21]. In particular, we were able to obtain
results on extremal problems related to the spaces V #,p (R d ) of functions of "bounded
variation" and Besov spaces B # (R d ). This algorithm is ready to implement for some
settings, depending on the value of p (if L p norm is chosen) and order of local polynomials
(it is more di#cult for r > 1), though the bookkeeping may be messy. On
the other hand, this algorithm is designed for the multivariate case. Its univariate
version would not only be much more complex than necessary, but would also produce
one-dimensional rings, that is, unions of the subintervals not necessarily neighboring,
which are unpleasant and, as it turned out, unnecessary. Our modified algorithms
take advantage of the simplicity of the real line topology, simply merging neighboring
intervals, thus resulting in partitions consisting of only intervals. These algorithms
cannot be easily generalized to multivariate setting, since a procedure of emerging
neighboring cubes may generate very complicated and undesirable sets in a partition.
This also makes it much more di#cult to establish Jackson inequalities for local ap-
proximants. We refer the interested reader to [21], where one can find that the proof
of Jackson inequality on a ring is already di#cult enough. For these reasons, we
strongly believe that simpler and more e#cient univariate algorithms are necessary.
2. Preliminaries. Throughout this paper, when we say that f, the function to
be approximated, belongs to L p (I), we mean f # L p (I) if 0 < p < #, and f # C(I)
r is an integer, 0 < # < r and 0 < p, q #, then the Besov space
is the set of all functions f # L p (I) such that the semi-(quasi)norm
sup
is finite, where # r is the usual rth modulus of smoothness. The (quasi)norm for
defined by
We also define a short notation for a special case that is used frequently in the theory:
If there is no potential confusion, especially in the case I = [0, 1], the interval I will be
omitted in the notation for the sake of simplicity. For example, L p stands for L p [0, 1]
are quasi-normed, complete, linear spaces continuously embedded
in a Hausdor# space X, then the K-functional for all f #
defined as
K(f, t, X 0 ,
This can be generalized if we replace # X1 by a quasi-seminorm | - | X1 on
K(f, t, X 0 ,
MODIFIED ADAPTIVE ALGORITHMS 1015
The interpolation space (X 0 , consists of all functions
< #, where
|f | (X0,X1 ) #,q
sup
When studying an approximation method, it is very revealing to know its approximation
spaces, which we now define. Let functions in a quasi-normed linear space X
be approximated by elements of its subsets # n , . , which are not necessarily
linear but are required to satisfy the assumptions
any a #= 0;
does not depend on n;
(v) # n=0 # n is dense in X;
(vi) Any f # X has a best approximation from each # n .
All approximant sets in this paper satisfy these assumptions. Denoting
we define the approximation space
A #
to be the set of all f # X for which E n (f) is of order n -# in the sense that the
following seminorm is finite:
|f | A #
sup
The general theorem below enables one to characterize an approximation space by
merely proving the corresponding Jackson and Bernstein inequalities (see [13, sections
7.5 and 7.9], [9], and [15]).
Theorem A. Let Y := Y # > 0, be a linear space with a semi-(quasi)norm | - | Y
that is continuously embedded in X. If {# n } satisfies the six assumptions above, and
Y satisfies the Jackson inequality
and the Bernstein inequality
then for all 0 < #, 0 < q #, the approximation space
A #
By a partition
of the interval [0, 1] we mean a finite set of subintervals
whose union # n
I . The (nonlinear)
Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
spaces # n,r of all PP functions of order r on [0, 1] with no more than n > 0 pieces
are defined by
I#P
I (x)# I (x), |P| # n},
where P I are in P r-1 , the space of polynomials of degree < r, and # I are the characteristic
functions on I. # 0,r is defined as {0}. These #
all assumptions (i)-(vi) on {# n } (see p. 3). The degree of best approximation of a
function f by the elements of # n,r is denoted by # n,r (f) p := E(f, # n,r ) p .
Remark . Some authors use the notation # (n-1)r,r in place of # n,r , since PP
functions can be viewed as special kinds of splines with each interior break point x i ,
a knot of multiplicity r. Also in use is PP n,r . Following
general notation in nonlinear approximation, we use the first subscript for the number
of coe#cients in the approximant. See [13], [14], [17], [26]. Strictly speaking, all n-
piece PP function of order r only form a proper subset of the free knot spline space
(n-1)r,r , but this subset has the same approximation power in L p as the whole space
(see Theorem 12.4.2 of [13]).
In his 1988 paper [23] (also see [24] and [13, section 12.8]), Petrushev characterized
the approximation space A #
using the Besov spaces; see the
following theorem.
Theorem B. Let 0 < p < #, n > 0, and 0 < # < r. Then we have
and
Therefore for 0 < q # and 0 < # < r
A #
In particular, if #
A #
The inequality (2.4) can be proved by finding a balanced partition
according to the function
s (f, x)| # dsdt
in the sense that
(see [13] for details of the proof). In fact, many Jackson-type inequalities can be
proved by showing the existence of a balanced partition (see, e.g., Theorems 12.4.3,
5, and 6 in [13], Theorem 1.1 in [19], and parts of Theorems 2.1 and 4.1 in [17]). We
state here Theorem 12.4.6 of [13], given by Burchard [8] in 1974 for the case
(see also de Boor [3]).
MODIFIED ADAPTIVE ALGORITHMS 1017
Theorem C. Let r and n be positive integers, and let # := (r
monotone function, then
#,p be the
space of functions f # L p [0, 1] for which the variation
|f | V #,p
I#P
is finite, where the sup is taken over all finite partitions P of [0, 1]. Following [17]
(see also Brudnyi [7] and Bergh and Peetre [1]), we define a modulus of smoothness
f,
0<h#t
sup
The following theorem, which is due to DeVore and Yu [17], provides characterization
of A #
using interpolation spaces involving V #,p .
Theorem D. Let 0 < p #, 0 < # < r, and #
approximation by elements from {# n,r } # 0 , we have the Jackson inequality
and the Bernstein inequality
Therefore
A #
In particular, if p < #,
A #
The Jackson inequality (2.12) follows from the definition of # f, t) #,p and the existence
for any f # V #,p of an S # n,r with n :=
which can be proved (see [17]) by showing the existence of a balanced partition
such that
and then defining S by
where P i are best L p approximations to f on I i from the space P r-1 .
Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
3. Adaptive algorithms.
3.1. The original adaptive algorithm. More than likely it will be hard to
find an exactly balanced partition numerically. An algorithm of this sort by Hu
[20], for instance, uses two nested loops (there is another level of loop that increases
the number of knots). This is probably one of the reasons why much attention is
paid to adaptive approximation, which selects break points by repeatedly cutting the
intervals into two equal halves, and produces PP functions with dyadic break points,
which can be represented by finite binary numbers of the form m - 2 -k ,
. Denote the spaces of such PP functions by # d
n,r and their approximation
errors E(f, # d
n,r (f) p . We now describe the original adaptive algorithms in
the univariate setting.
Let E be a nonnegative set function defined on all subintervals of [0, 1] which
satisfies
E(I) # E(J) if I # J ;
uniformly as |I| # 0.
Given a prescribed tolerance # > 0, we say that an interval I is good if E(I) #;
otherwise it is called bad. We want to generate a partition G := G(#, E) of [0, 1] into
good intervals. If [0, 1] is good, then is the desired partition; otherwise
we put [0, 1] in B, which is a temporary pool of bad intervals. We then proceed with
this B and divide every interval in it into two equal pieces and test whether they are
good, in which case they are moved into G, or bad, in which case they are kept in B.
The procedure terminates when resulting intervals are good
and are in G), which is guaranteed to happen by (3.2).
The set function E(I) usually depends on the function f that is being approximated
and measures the error of approximation of f on I, such as # I G(x) dx in (2.9),
thus will be called the (error) measure of I throughout this paper. In the simplest
case, E(I) is taken as the local approximation error of f on I # [0, 1] by polynomials
of degree < r:
E(I)
and the corresponding approximant on G is defined by (2.15). This gives an error
where |G| is the number of intervals in G. One can estimate in di#erent ways
a n (f) p := a n (f, E) p := inf |G| 1/p #,
where the infimum is taken over all # > 0 such that
and Solomjak [2] and DeVore [12] for estimates for functions f in Sobolev spaces.
Other estimates can be found in Rice [25], de Boor and Rice [6], and DeVore and
Yu [18] and the references therein. We only mention the following two results.
Theorem E (see [18, Theorem 5.1]). Let
If f # C r (0, 1) with |f (r) (x)| #(x), where # L # is a monotone function such that
I
MODIFIED ADAPTIVE ALGORITHMS 1019
where C 1 is an absolute constant, then we have
a n (f) # Cn -r # .
Note that compared with Theorem C with
theorem has an extra requirement (3.5) on #.
Theorem F (see [18, Corollary 3.3]). Let 0 < p < # > 0, and q > # :=
a
# (L q ) , we see (3.7) is weaker than (2.4), which is for free
knot spline approximation. The reason for this is not hard to see: adaptive algorithms
not only select break points from a smaller set of numbers (that is, the set of all finite
binary numbers), but they also do it in a special order. Consider
as an example, a good free knot approximant will have most knots very close to 0 (see
examples in [20] and Table 5.2 later in this paper). However, an adaptive algorithm
needs at least n - 1 knots, 2 before it can put one at 2 -n and thus
needs more knots than a free knot spline algorithm. Although one classifies adaptive
approximation as a special kind of free knot spline approximation (since the knots
sequence depends on the function to be approximated), one is far from free when
choosing knots. It is considered "more restrictive" (DeVore and Popov [14]) than free
knot spline approximation.
We should point out that all theorems mentioned in this subsection are of a
Jackson-type, that is, so-called direct theorems. Bernstein inequalities (closely related
to inverse theorems, sometimes referred to also as inverse theorems themselves) for free
knot splines, such as (2.5) and (2.13), are valid for all splines, including PP functions
produced by adaptive algorithms. The problem is that all Jackson inequalities for
the original adaptive algorithms are not strong enough to match those Bernstein
inequalities in the sense of Theorem A. From this point of view, Theorems E and
F are weaker than they look. We do not know exactly what kind of functions can
be approximated by the original adaptive algorithms to a prescribed order, that is,
we can not characterize their approximation spaces A #
q . They do not fully exploit the
power of nonlinear approximation, and sometimes they generate too many intervals,
many of which may have an error measure much smaller than #.
As mentioned above, there are two major aspects in which adaptive approximation
is di#erent from free knot spline approximation: (a) a smaller set of numbers to choose
knots from and (b) a special, and restrictive, way to select knots from the set. It turns
out that (b) is the reason for its drawback. Although it is also the reason why adaptive
approximation is simple (and we want to keep it that way), it does not mean we have to
keep all the knots it produces. In this paper, we modify the usual adaptive algorithm
in two ways. The idea is that of splitting AND merging intervals/cubes used in a
recent paper by Cohen et al. [11]. The two new algorithms generate partitions of
[0, 1] with fewer dyadic knots which are nearly balanced in some sense. In section 4,
we prove that they have the same approximation order as that of free knot splines.
3.2. Algorithm I. We start with the original adaptive procedure with some
# > 0, which generates a partition
of [0, 1] into good intervals. The
number N # may be much larger than it has to be. To decrease it, we merge some
of the intervals I #
i . We begin with I #
1 and check the union of I #
1 and I # 2 . If it is still
a good interval, that is, if its measure E(I #
#, we add I # 3 to the union and
1020 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
check whether E(I #
3 ) #, and we proceed until we find the largest good union
I #
k in the sense that
but
E(I #
We name I # 1 # I #
k as I 1 . If k < N # , we continue with I #
k+1 and find the next
largest good union as I 2 . At the end of this procedure, we obtain a modified partition
consisting of N # N # good intervals
for which each union J i := I i # I i+1
is bad,
This partition is considered nearly balanced. For the size of N we have
# .
3.3. Algorithm II. Our second algorithm generates a nearly balanced partition
in another way. It does not make heavy use of prescribed tolerance #; rather, it merges
intervals with relatively small measures while dividing those with large ones. As in
the ordinary adaptive algorithms, we start with dividing [0, 1] into two intervals I 1
and I 2 of equal length. However, this is where the similarity ends. We then compare
measures E(I 1 ) and E(I 2 ) and divide the interval with larger measure into two equal
pieces. In the case of equal measure, we divide, rather randomly, the one on the left.
Now we have three intervals and are ready for the three-step loop below.
Step 1. Assume there is currently a partition {I i } k
I j has the largest
measure among all I i . If E(I j+1
is a fixed parameter, we check the union of I j+1 # I j+2 to see whether its measure
E(I j+1 # I j+2 ) < M . If so, add the next interval I j+3 into the union and check its
measure again. We continue until we get a largest union # j+m1
I i whose measure is
less than M, and replace this union by the intervals it contains. Then, if j +m 1 < k,
we find the next largest union # j+m1+m2
I i in the same manner and replace these
intervals by their union. Furthermore, we do the same to the intervals to the left of
I j (but keep I j intact). In this way we obtain a new partition with (the old) I j still
having the largest measure. This partition is nearly balanced in the sense that the
measure of the union of any two consecutive new intervals is no less than #
(because these new intervals were largest unions of old intervals). At the end of this
step we renumber the new intervals and update the value of k.
Step 2. Check whether the new partition produced in Step 1 is satisfactory using
an application-specific criterion, for instance, whether k has reached a prescribed value
n or the error is reduced to a certain level. If not, continue with Step 3; otherwise
define the final spline by (2.15) and terminate the algorithm.
Step 3. Divide the interval with the largest measure into two equal pieces,
renumber the intervals, update the values of k and M, and then go back to Step 1.
Remark. In Step 1, if I l and I l+1 are the two newest intervals (two "brothers"
with equal length), one needs only to check I l-1 # I l if l - 1, l #= j, and/or I l+1 # I l+2
since other unions of two consecutive intervals have measures no
MODIFIED ADAPTIVE ALGORITHMS 1021
less than the value of M in the previous iteration, which is, in turn, no less than the
current M . We stated it in the way above only because it shows the purpose of the
step more clearly.
It should be pointed out that one needs to be careful about the stopping criterion
in Algorithm II. For example, if it is applied to the characteristic function
after two iterations we will always
have The break point # 2/2 in this example can be
replaced by any number in (0, 1) which does not have a finite binary representation
such as 0.4. If k # used as the sole stopping criterion, the algorithm will fall
into infinite loop. Fortunately, the error in this example still tends to 0; therefore,
infinite loop can be avoided by adding error checking in the criterion. The next lemma
shows this is the case in general.
Lemma 3.1. Let E be an interval function satisfying (3.1) and (3.2), and let
prescribed. Then the criterion
will terminate Algorithm II.
Proof. We show that if k never exceeds n, then as the number
of iterations goes to #. Let 0 < # < 1 be fixed. Let
with the
max taken at one moment. Fix this
M and denote the group of all subintervals in
the partition with "large" errors by G
be as in Step 1, changing from iteration to iteration. We have
# M from now on.
We first make a few observations. Since the interval currently having the largest
measure is always in G
M , each iteration cuts a member of G
M . However, the algorithm
will not merge any member I i # G
M with another interval because E(I i
any union of I i with another interval would have even larger measure by
(3.1). By (3.2), there exists # > 0 such that |I i | > # for any I i # G
M . Note all
intervals in a partition are disjoint, thus the total length of the intervals in G
M is no
larger than 1, and its cardinal number |G #
M | # 1/#.
From these observations, we conclude the following. When an iteration cuts a
member I i of G
M into two "children" of equal length, one of the three cases will
happen: (a) neither child of I i belongs to G
M , thus |I i | > # is removed from the total
length of G
exactly one of the children belongs to it (hence having a length
> #) and the other child, with the same length |I i |/2 > #, is removed from G
or (c)
both children belong to it. The case (a) decreases |G #
M | by 1, (b) keeps it unchanged,
and (c) increases it by 1. Now one can see that at most # 3/# +1 iterations will empty
G
M , since at least one third of them will be cases (a) or (b) to keep |G #
M | # 1/#, which
will remove all the total length of G
M , thus emptying it. This reduces the maximum
error by a factor # < 1. Repeat this enough times and the maximum error
will eventually tend to 0.
Although (3.2) does not say anything about the convergence rate of E(I) as |I| #
0, and the proof of the above lemma may make it sound extremely slow, one can expect
a fairly fast convergence in most cases. For example, in the case
if f is in the generalized Lipschitz space Lip #, p) := Lip #, L p [a, b]), 0 < # < r,
that is, if
|f | Lip # := |f | Lip #,p) := sup
1022 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
then for any I # [a, b]
|f | Lip # .
We feel it is safe to say that most functions in applications belong to Lip #, p) with
an # reasonably away from 0, at least on subintervals not containing singularities,
thus halving an interval often reduces its error by a factor of 2 # .
A natural question that may arise here is: How complex are the new algorithms?
We give brief comparisons below to answer this question. Algorithm I is straight for-
ward. It is the original adaptive algorithm with a second (merging) phase added.
This phase consists of no more than merging attempts, where N # is
the number of subintervals the original algorithm generates, and N that of the final
subintervals. As for Algorithm II, there are two major di#erences from the original
version. The first one, as mentioned in the remark after the algorithm description,
is: up to two merging attempts are made after cutting each interval. The other one
is in the book-keeping. In the original version, a vector is needed to record errors
on all intervals (or to indicate which intervals are bad), while Algorithm II keeps the
index of the interval that has the largest error E(I) in a scalar variable, in addition
to the vector containing all errors. This requires a search for the largest element in
the vector after each cutting or merging operation.
One can see from above that the new algorithms are not much more complex
in terms of programming steps. The added CPU time, in terms of the number of
results mainly from the evaluations of the error measure E(I)
required by merging operations. Our estimate is that either algorithm uses two or
three times as much CPU time as the original algorithm. More information on CPU
time will be given in section 5 together with numerical details.
4. Approximation power of the algorithms. We now show that our modified
adaptive algorithms have the full power of nonlinear approximation. More precisely,
we prove that they produce piecewise polynomials satisfying the very same Jackson
inequalities for free knot spline approximation (with possibly larger constants on the
right-hand side since the partitions are not exactly balanced). As we mentioned
earlier, the corresponding Bernstein inequalities hold true for all splines; therefore we
are really proving that the approximation spaces for the modified adaptive algorithms
are the same as those for free knot spline approximation.
We state below our results as three main theorems, parallel to Theorems B, C,
and D, respectively. In fact, we can prove most results of this kind for our algorithms,
such as Kahane's theorems and its generalization [13, Theorems 12.4.3 and 5], but
the proofs would be too similar to the ones below.
We recall that throughout this paper, I j denotes the interval with largest measure
among all I i in the partition, the union of any two consecutive intervals J
has a measure E(J i ) > E(I j ), and J i is called bad in Algorithm I. All PP functions
on the resulting partitions are defined by (2.15).
Theorem 4.1. Let n and r be positive integers, and let 0 < p < #, 0 < # < r,
then the two modified adaptive algorithms (with
defined in (2.8) or (ii)
functions S of (2.15) that satisfy
the Jackson inequality
(4.
MODIFIED ADAPTIVE ALGORITHMS 1023
From Theorem A we obtain the approximation space A #
product. It turns out to be the same as A #
which is not surprising
since # d
n,r is dense in # n,r . The surprising part is that one can get such an approximant
using a simple adaptive algorithm.
Corollary 4.2. Let 0 < p < #, 0 < q #, 0 < # < r, and
. For approximation by PP functions in # d
n,r , we have
A #
In particular,
A #
Proof of Theorem 4.1. The proofs of the theorem in the cases (i) and (ii) are
very similar. We only consider (i) and remark that, in the case (ii), the inequality
plays the major role.
PP approximants produced by Algorithm I. Let E(I) := # I G(x) dx, where G is as
in (2.8), and # :=
We claim that the number N of intervals it produces is no greater than 2n+1. Indeed,
by (3.8)
The rest of the proof of (4.1) is similar to that of (2.4) (cf. section 12.8, p. 386 of
[13]); we sketch it here for completeness. It is proved in [13] that for any f # B # [0, 1],
M is equivalent to |f | B # [0, 1] with constants of equivalence depending only on r and
#, and that for such an f
Define the approximant S by (2.15) and we have
here in the fifth step we have used the equality # and in the last step
we have used the equivalence of M and |f | B # .
PP approximants produced by Algorithm II. Let E(I), M, and # be the same as
above, and use (3.9) as stopping criterion in Step 2. If the algorithm terminates due
1024 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
to (thus giving less than n pieces), it is the same situation as with
Algorithm I. Otherwise we have n pieces when it terminates, and (4.1) follows:
Cn
Theorem 4.3. Under the conditions of Theorem C, the modified adaptive algorithms
(with E(I) := # I #(x) # dx and # := n
produce PP approximants
S in # d
n,r that satisfy the Jackson inequality:
Proof of Theorem 4.3.
PP approximants produced by Algorithm I. Let
E(I) := # I
as the Taylor polynomial for f of degree r - 1 at the point x i+1 (not best
we have (see equation (4.15) in Chapter 12 of [13])
Using (4.6) in place of (4.4), then (4.5) for p < # can be proved by arguments very
similar to those in the proof of Theorem 4.1 by Algorithm I. We also refer the reader
to the proof of Theorem C in [13]. For #, the estimate of N is the same and we
need only to replace # N
by
PP approximants produced by Algorithm II. Let E(I), #, and M be the same as
above. Use (3.9) again as the stopping criterion in Step 2. If the algorithm terminates
because it is the same situation as in Algorithm I. Otherwise, for
p/#
MODIFIED ADAPTIVE ALGORITHMS 1025
where we have used the inequality (4.6) in the second step, and
the last one. For #, we make similar changes to those in Algorithm I:
Theorem 4.4. Let n and r be positive integers, and let 0 < p, #, 0 < # < r,
then the two modified adaptive algorithms (with
E(I)
#,p ) produce PP functions S of (2.15) that
satisfy the Jackson inequality
Using Theorems 4.4 and A we have the following characterization of A #
Corollary 4.5. For approximation by PP functions in # d
n,r we have
A #
In particular, if p < #,
A #
Proof of Theorem 4.4. It su#ces to show (2.14) since (4.7) immediately follows
from it with any t > 0 and n := (see the end of section 2). We only prove it
for p < #. The case of can be verified by making changes similar to those in
the proof of the L# case in the previous theorem.
PP approximants produced by Algorithm I. Let E(I)
p and
#,p . From (3.8), the number N of intervals the algorithm produces
can be estimated as
Indeed, if N > 2n (otherwise, it's done) we have
#,p
Cn #/p
#,p
where we have used the definition (2.11) of # f, t) #,p . Since 1 - #/p, this gives
Cn. Now (2.14) follows, since
#,p .
PP approximants produced by Algorithm II. We set E(I) := E r (f, I) #
, and use
#,p in the stopping criterion (3.9). If it stops
because have exactly the same situation as with Algorithm I, with
1026 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
the same partition; otherwise there are n intervals when it terminates. In the latter
case, we have
5. Numerical implementation and examples. Theoretically, the two algorithms
have the same approximation power. However, when it comes to numerical
implementation, we prefer Algorithm II since it directly controls the number of polynomial
pieces n, while # in Algorithm I is neither a power of n nor a tolerance for
(though it is closely related to both). We implemented Algorithm II on the
computer, using Fortran 90 and mainly for 2. The error measure used in the code
is
2 unless we have a better one to use, such as # I # I |f (r)
| # for
the square root function in the first example in this section. The L 2 norm of f on the
interval I estimated by the composite Simpson rule for integral
and its L# norm is estimated by
are equally distributed nodes, n p is a program
parameter roughly set as 6 times r, and . The best L 2 polynomial
approximant on I i , discretized by (5.1) as an overdetermined n p - r system of linear
equations for the least squares method, is calculated by either QR decomposition or
singular value decomposition by calling LINPACK subroutines Sqrdc and Sqrsl, or
Ssvdc (or their double precision counterparts). The latter takes longer but we did
not see any di#erence in the first four or five digits of the local approximation errors
they computed; thus we did not test it extensively.
The L# version of algorithm is basically the same, except that we use
estimated by (5.2). The local polynomials P I (and the global smooth splines)
are still obtained by the least squares method, that is, still best L 2 approximants.
This is common in the literature, and it is justified by the fact that the best L 2
polynomial approximant is also a near-best L# polynomial on the same interval; see
Lemma 3.2 of DeVore and Popov [16].
The number of polynomial pieces is used as the main termination criterion, while
# in (3.9) is set to a small value mainly to protect the program from falling into infinite
loops, rather than the sophisticated ones as in proofs in the previous section. It turned
out that infinite loop is not a problem. A nonfull rank matrix in the least squares
method is a problem, which happens far before it falls into an infinite loop. This
is because if I i is too small, the machine will have di#culties distinguishing the n p
MODIFIED ADAPTIVE ALGORITHMS 1027
points needed in (5.1). Therefore, we added a third condition to protect the program
from failing: stop the program when
We also added a second part in the code, namely, finding an L 2 smooth spline
approximation to the function with the knot sequence {t i } n+r
i=2-r , where the interior
knots a < t 2 < t 3 < - < t n < b are the break points of the PP function obtained
by Algorithm II, used as single knots, and the auxiliary knots are set as t
b. Despite the fact that the
partitions are guaranteed to be good only for PP functions, they usually work well
for smooth splines, too. De Boor gave some theoretical justification in the discussion
of his subroutine Newnot [4, Chapter XII].
The least square objective function for finding this smooth spline -
S is
set as 5r+1, is the number of equal pieces into which we cut each subinterval
I are the points resulted from such cutting, and the weights w j are chosen so
that (5.4) becomes a composite trapezoidal rule for the integral # b
a # f(x) -
dx:
The actual calculation
of the B-spline coe#cients of
are the B-splines with the knot sequence {t i } scaled
so that
done by de Boor's subroutine L2Appr in [4, Chapter XIV].
We used the source code of all the subroutines in the book from the package PPPACK
on the Internet.
We tested our code on a SUN UltraSparc, with a clock frequency 167MHz, 128MB
of RAM, and running Solaris 2.5.1. The speed is so fast that it is not an issue here:
for finding break points, it is somewhere from 0.015 second for to 0.1 second for
printing minimum amount of messages on the screen, and it is less than
10% of these for computing smooth splines. We also tested the code on a 300 MHz
Pentium II machine with 64 MB of RAM running Windows NT 4.0. The speed is at
least three times as fast. None of the problems we tested used more than 0.1 second.
(The reason for the great di#erence in speed may be that the SUN we used is a file
server, not ideal for numerical computation.) There is still room for improvement in
e#ciency. For example, one can use a value of n p , larger than what we use, at the
beginning and decrease it as n increases (and the error on each subinterval decreases).
The value of n s should be related to n, too, for the same reason.
The main cost of CPU time is the evaluation of the error measure E(I) for each
subinterval I. We use
estimated by QR decomposition, as an exam-
ple. Each such problem involves n p function evaluations, and (n p - r
operations required in QR decomposition, plus some more for estimating the error
from the resulting matrices. Each cutting of intervals requires two E(I) evaluations,
1028 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU

Table
Approximation order of
and each merging attempt requires one. Our numerical experiments show that a typical
run resulting in n subintervals cuts intervals about 2n time. Each cutting results
in up to two attempts of merging subintervals. That gives about 8n least squares
problems, each of which involves n p function evaluations plus about n p r 2 arithmetic
operations. In view of the approximation order we proved in the previous section, and
the fact that n p is roughly a multiple of r, we think it pays to use a relatively large r,
at least 4 or 5. For 5, the error will reach the machine epsilon (single precision)
when n is somewhere between 30 and 70 in most cases.
We use the square root function to test the PP function approximation
order. This function is only in the Lipschitz space Lip( 1
thus the approximation
order is only 1/2 for splines with equally spaced knots in the L# norm, no matter
what their order r is. By Theorem 4.3, we should have e n := #f -S n # p # Cn -r , where
S n is the function consisting of n polynomial pieces computed by Algorithm II using
I |f (r) (x)| # dx, and we have combined # in the theorem
into the constant C. After the knot sequence has been found, QR decomposition is
used at the end of the program on each subinterval to estimate e n . Since the error
decreases fast for double precision had to be used in QR decomposition for
large values of n. Assume that what we actually obtain from the code is e
where # is the approximation order. Since log e plot the
points in the plane, they should form a line. Since such
a plot zigzags very much, we calculated the least squares line fitting
to find the order. Table 5.1 gives values of # for di#erent r using both L 2 and
L# norms. We should mention that the points values of n are too
low and ruin the obvious line pattern formed by those for larger n, thus we give two
values of #, one from the points for and the other from
As can be seen from the table, the latter values are right around or even exceed r.
Remark. We tried some power of E r (f, I) p for E(I) and felt, in view of (4.6),
it would yield a better balance of subintervals, thus a higher order. But the orders so
obtained were well below r (4.46 for e.g. The reason might be that
I # is additive, but (power of) E r (f, I) p is not.
To illustrate the advantage of interval merging, we compare the original adaptive
algorithm and our modified ones with the function
log 2
-m
This function is in C # , and is decreasing and convex on [0, 1] with
. Note that
since f is decreasing on [0, 1]. Table 5.2 shows
comparison in numbers of knots produced for the same approximation error by the
original adaptive algorithm and our Algorithm II. Both programs try to put first knots
near where the graph is very steep. The original algorithm has to, as pointed
out early, lay down knots 2 one by one before reaching an error of
MODIFIED ADAPTIVE ALGORITHMS 1029

Table
Comparison in numbers of interior knots produced by the original and modified adaptive algorithms
for the same error in approximating
Original
Alg.
while Algorithm II, after trying all these knots one at a time and merging all
but the last interval, puts the very first knot at 2 -23 .
It is interesting to watch how Algorithm II moves a knot toward a better position
in successive iterations without increasing the total number of pieces. The following
screen output shows that in iterations 1 and 2 the program moves the break point 0.5
to 0.25 and then to 0.125, while the error decreases form 0.5 to 0.47; in iterations 3-22
it moves the break point all the way to 2 -23 # with the error decreased
to 0.27. What happened internally is, in iteration 1, e.g., it cuts the interval [0, 0.5]
into [0, 0.25] and [0.25, 0.5]. Since the error on the union of [0.25, 0.5] and [0.5, 1] is
smaller than that on [0, 0.25], it then merges the two intervals into [0.25, 1]. The net
e#ect of these steps is moving the break point 2 -1 to 2 -2 .
Iteration 0: # of
errors=
L_\infty error on [a,
Iteration 1: # of
errors=
L_\infty error on [a,
Iteration 2: # of
errors=
L_\infty error on [a,
(Many lines deleted.)
Iteration 22: # of
errors=
2.70000E-01 2.30000E-01
L_\infty error on [a,
Y.-K. HU, K. A. KOPOTUN, AND X. M. YU

Table
Approximation errors to the Runge function on [-5, 5].
9
We now consider the infamous Runge function, which is also in C # but, on the
other hand, is hard to interpolate or approximate. Lyche and M-rken [22] approximated
it by the knot removal algorithm, and Hu [20] approximated it by balancing
the rth derivative of the function on subintervals in two nested loops. Here and in
the rest of the paper, we use 4. In Table 5.3, we compare our results with those
of Lyche and M-rken (LM) [22] and Hu [20]. For the same number of knots (that is,
list our errors measured in # 2 / # b - a for the PP function S n and the
smooth spline -
also that of -
measured in L# norm. We divide the L 2 norm by
since it is more comparable to the L# norm, which is what LM and Hu used.
The errors by LM are estimated from figures in [22]. Because of the simple nature of
our algorithm, we only expected to compete with their results by splines with two or
three times as many knots. It turns out that our approximation errors are almost as
good as theirs, which were produced by more sophisticated methods.
By now, the reader may begin to wonder: what is the e#ect of the parameter
of Algorithm II, used in Lemma 3.9 to guarantee the termination
of Algorithm II. We tried functions we tested, it worked excellently
except that the number of polynomial pieces went up and down a few times with the
square root function using dx, in which case used
instead. It is true that in theory it might get into an infinite loop, but since our goal
is to find a nearly balanced partition, better in this aspect, provided
infinite loop does not happen. It did not. As a matter of fact, sometimes we feel
the need for a value slightly larger than 1, e.g., with symmetric functions such as the
Runge function. What happens with # 1 is that if there are two subintervals having
the same largest measure at the moment, symmetric about the center of the interval,
then the outcome of the next iteration, which processes the subinterval on the left,
will very often interfere with the processing of the subinterval on the right later. It
may not make the approximation error worse, at least not by much, it is just that the
knot sequence becomes unsymmetrical, thus unnatural and unpleasant. Furthermore,
most algorithms in the literature produce symmetric knots for symmetric functions;
it would be hard to compare our results with theirs. For these minor reasons, we
used in preparation of Table 5.3. In the next example, we consider the PP
function
which has a jump at # 2/2. As we mentioned in the discussion before Lemma 3.1, since
# 2/2 has no finite binary representation, this function can never be approximated
exactly by a PP function with dyadic break points. The program (with
cutting and merging around the jump (since the number of pieces is always 3 after
two iterations), until it is stopped by the criterion (5.3), resulting in t
MODIFIED ADAPTIVE ALGORITHMS 1031
500 600 700 800 900 1000 11000.611.41.82.2-.04
Temperature
Fig. 5.1. Titanium Heat Data (circles). The final spline (solid line) has 15 interior knots. The
errors for preapproximation (dotted) and for the final spline (dashed) use scales on the right.
and 0.70710754. The PP function matches f exactly on the computer screen
since the two points are indistinguishable. One can very well combine them into a
single break point, thus virtually reproducing f . The original adaptive algorithm, in
contrast, would put many many knots around the jump while trying to narrow the
subinterval containing the jump: 0.5, 0.75, 0.625, 0.6875, 0.71875, . All these knots
are useless except the newest two.
In practice, one often wants to approximate discrete data points other than known
functions as in the previous examples. In this case, we preapproximate the points by
a spline with as many parameters as we wish to use, then apply our algorithm to this
spline. For smooth-looking data, we interpolate the data by a C 1 cubic spline with
knots at the data points, using de Boor's subroutine Cubspl in [4]. This worked very
well. We produced some sample data points from the Runge function and square root
function and applied this approach to them. It resulted in virtually the same knot
sequences as those generated by directly approximating the original functions.
In the real world, however, it is likely that the data will contain errors. If the
data points are interpolated, one can see small wiggles in the graph, which tricks the
program laying knots in areas where the curve is otherwise flat. One such example
is the Titanium Heat Data (experimentally determined), see [4, Chapter XIII], and
also LM [22] and Hu [20]. In Figure 5.1 the reader can see wiggles on both the left
and right. De Boor [4, Chapter XIV] suggests that the data be approximated by a
less smooth spline. We absolutely agree. For the same reason, we used fewer knots
for preapproximating spline in the flat parts at both ends, than we did near the high
1032 Y.-K. HU, K. A. KOPOTUN, AND X. M. YU
peak around 900 # , trying to ignore the wiggles. In fact, we used almost the same knot
sequence for preapproximating spline as in Figure 4 of [20].

Table
Approximation errors to the Titanium Heat Data.
Obtained by Order # of knots Error
Alg. II 4 11 0.070
Alg. II 4 15 0.031
Since de Boor, LM, and Hu all used L# norm for approximating these data, we
also used the L# version of our program. Figure 5.1 shows a cubic spline approximation
to the Titanium Data obtained by this method. It has 15 interior knots with
an error of 0.031. Table 5.4 gives a comparison of our results with those by others on
the same data.

Acknowledgments

. We are deeply indebted to Professor Ron DeVore, who
inspired us by discussing the excellent ideas in [11] during our visit to the University
of South Carolina. We want to thank him and Professors Pencho Petrushev and
Albert Cohen for providing us with drafts of their manuscript [11]. Credit is also due
to Professor Dietrich Braess, the editor of this paper, and the referees, whose opinions
and suggestions helped very much in improving the manuscript. As a matter of fact,
we reshaped the last section during the communication with them.



--R

On the space Vp (0
Piecewise polynomial approximation of functions of classes W
Good approximation by splines with variable knots
A Practical Guide to Splines
Least squares cubic spline approximation II-Variable knots
An adaptive algorithm for multivariate approximation giving optimal convergence rates
Spline approximation and functions of bounded variation
Splines with optimal knots are better
Jackson and Bernstein-type inequalities for families of commutative operators in Banach spaces
Adaptive wavelet methods for elliptic operator equations-Convergence rates
Nonlinear approximation and the space BV (R 2
A note on adaptive approximation


in Function Spaces and Applications
Interpolation of Besov spaces

Degree of adaptive approximation
Convexity preserving approximation by free knot splines
An algorithm for data reduction using splines with free knots
On multivariate adaptive approximation
A data reduction strategy for splines with applications to the approximation of functions and data
Direct and converse theorems for spline and rational approximation and Besov spaces
Rational Approximation of Real Functions

Basic Theory
--TR
adaptive algorithms;nonlinear approximation;data reduction;besov spaces;degree of approximation;piecewise polynomials;splines;modulus of smoothness;approximation spaces