Learning and Design of Principal Curves.
AbstractPrincipal curves have been defined as self-consistent smooth curves which pass through the middle of a d-dimensional probability distribution or data cloud. They give a summary of the data and also serve as an efficient feature extraction tool. We take a new approach by defining principal curves as continuous curves of a given length which minimize the expected squared distance between the curve and points of the space randomly chosen according to a given distribution. The new definition makes it possible to theoretically analyze principal curve learning from training data and it also leads to a new practical construction. Our theoretical learning scheme chooses a curve from a class of polygonal lines with $k$ segments and with a given total length to minimize the average squared distance over $n$ training points drawn independently. Convergence properties of this learning scheme are analyzed and a practical version of this theoretical algorithm is implemented. In each iteration of the algorithm, a new vertex is added to the polygonal line and the positions of the vertices are updated so that they minimize a penalized squared distance criterion. Simulation results demonstrate that the new algorithm compares favorably with previous methods, both in terms of performance and computational complexity, and is more robust to varying data models.
Introduction
Principal component analysis is perhaps the best-known technique in multivariate analysis
and is used in dimension reduction, feature extraction, and in image coding and enhancement.
Consider a d-dimensional random vector moments.
The first principal component line for X is a straight line which has the property that the
expected value of the squared Euclidean distance from X to the first principal component
line is minimum among all straight lines. This property makes the first principal component
a concise one-dimensional approximation to the distribution of X, and the projection of
X to this line gives the best linear summary of the data. For elliptical distributions the
first principal component is also self consistent, i.e., any point of the line is the conditional
expectation of X over those points of the space which project to this point.
Hastie [1] and Hastie and Stuetzle [2] (hereafter HS) generalized the self consistency
property of principal components and introduced the notion of principal curves. Let
be a smooth (infinitely differentiable) curve in R d parametrized by t 2 R,
and for any x 2 R d let t f
(x) denote the parameter value t for which the distance between x
and f(t) is minimized (see Figure 1). More formally, the projection index t f (x) is defined by
denotes the Euclidean norm in R d .
Figure

1: Projecting points to a curve.
By the HS definition, the smooth curve f(t) is a principal curve if
f does not intersect itself
(ii) f has finite length inside any finite ball of R d
(iii) f is self-consistent, i.e.,
Intuitively speaking, self-consistency means that each point of f is the average (under the
distribution of X) of points that project there. Thus, principal curves are smooth self-consistent
curves which pass through the "middle" of the distribution and provide a good
one-dimensional nonlinear summary of the data.
Based on the self consistency property, HS developed an algorithm for constructing principal
curves. Similar in spirit to the Generalized Lloyd Algorithm (GLA) of vector quantizer
design [3], the HS algorithm iterates between a projection step and an expectation step.
When the probability density of X is known, the HS principal algorithm for constructing
principal curves is the following.
be the first principal component line for X. Set
\Psi for all x 2 R d .
Step 4 Evaluate the expected squared distance \Delta(f (j)
and f (j) . Stop if \Delta(f (j) ) is less than a certain threshold. Otherwise, increase j by 1
and go to Step 1.
In practice, the distribution of X is often unknown, but a data set consisting of n samples
of the underlying distribution is known instead. In the HS algorithm for data sets, the
expectation in Step 1 is replaced by a smoother (locally weighted running lines [4]) or a
nonparametric regression estimate (cubic smoothing splines). HS provide simulation examples
to illustrate the behavior of the algorithm, and describe an application in the Stanford
Linear Collider Project.
Alternative definitions and methods for estimating principal curves have been given subsequent
to Hastie and Stuetzle's groundbreaking work. Banfield and Raftery [5] modeled
the outlines of ice floes in satellite images by closed principal curves and they developed a
robust method which reduces the bias and variance in the estimation process. Their method
of clustering about principal curves led to a fully automatic method for identifying ice floes
and their outlines. On the theoretical side, Tibshirani [6] introduced a semiparametric model
for principal curves and proposed a method for estimating principal curves using the EM
algorithm. Close connections between principal curves and Kohonen's self-organizing maps
were pointed out by Mulier and Cherkassky [7]. Recently, Delicado [8] proposed yet another
definition based on a property of the first principal components of multivariate normal
distributions.
There remains an unsatisfactory aspect of the definition of principal curves in the original
HS paper as well as in subsequent works. Although principal curves have been defined to
be nonparametric, their existence for a given distribution or probability density is an open
question, except for very special cases such as elliptical distributions. This also makes it
very difficult to theoretically analyze any estimation scheme for principal curves.
In this paper we propose a new definition of principal curves to resolve this problem.
In the new definition, a principal curve is a continuous curve of a given length L which
minimizes the expected squared distance between X and the curve. In Section 2 (Lemma 1)
we prove that for any X with finite second moments there always exists a principal curve
in the new sense. We also discuss connections between the newly defined principal curves
and optimal vector quantizers. Then we propose a theoretical learning scheme in which the
model classes are polygonal lines with k-segments and with a given length, and the algorithm
chooses a curve from this class which minimizes the average squared distance over n training
points. In Theorem 1 we prove that with k suitably chosen as a function of n, the expected
squared distance of the curve trained on n data points converges to the expected squared
distance of the principal curve at a rate O(n 1=3 ) as n !1.
Two main features distinguish this learning scheme from the HS algorithm. First, the
polygonal line estimate of the principal curve is determined via minimizing a data dependent
criterion directly related to the definition of principal curves. This facilitates the theoretical
analysis of the performance. Second, the complexity of the resulting polygonal line is
determined by the number of segments k, which is typically much less than n for the optimal
choice of k. 1 This agrees with our mental image that principal curves should provide a
concise summary of the data. On the other hand, for data points the HS algorithm with
scatterplot smoothing produces polygonal lines with n segments.
Though amenable to analysis, our theoretical algorithm is computationally burdensome
for implementation. In Section 3 we develop a suboptimal algorithm for learning principal
curves. The practical algorithm produces polygonal line approximations to the principal
curve just as the theoretical method does, but global optimization is replaced by a less
complex iterative descent method. In Section 4 we give simulation results and compare our
algorithm with previous work. In general, on examples considered by HS, the performance
of the new algorithm is comparable with the HS algorithm, while it proves to be more robust
to changes in the data generating model.
We note here that the choice of k can be made automatic in principle by using the method of structural
risk minimization [9].
Learning Principal Curves with a Length Constraint
A curve in d-dimensional Euclidean space is a continuous function f : I ! R d , where I is
a closed interval of the real line. Let the expected squared distance between X and f be
defined by
\Theta inf
where the projection index t f
(x) is given in (1). Let f be a smooth (infinitely differentiable)
curve and for - 2 R consider the perturbation f + -g of f by a smooth curve g such that
proved that f is a principal curve if and only if f
is a critical point of the distance function in the sense that for all such g,
@-
0:
It is not hard to see that an analogous result holds for principal component lines if the
perturbation g is a straight line. In this sense the HS principal curve definition is a natural
generalization of principal components. Also, it is easy to check that principal components
are in fact principal curves if the distribution of X is elliptical.
An unfortunate property of the HS definition is that in general it is not known if principal
curves exists for a given source density. To resolve this problem we go back to the defining
property of the first principal component. A straight line s(t) is the first principal component
if and only if
\Theta min t
\Theta min t

for any other straight line - s(t). We wish to generalize this property of the first principal
component and define principal curves so that they minimize the expected squared distance
over a class of curves rather than only being critical points of the distance function. To do
this it is necessary to constrain the length 2 of the curve, since otherwise for any X with
a density and any ffl ? 0 there exists a smooth curve f such that \Delta(f ) - ffl, and thus a
minimizing f has infinite length. On the other hand, if the distribution of X is concentrated
on a polygonal line and is uniform there, the infimum of the distances \Delta(f ) is 0 over the
class of smooth curves, but no smooth curve can achieve this infimum. For this reason, we
2 For the definition of length for nondifferentiable curves see Appendix A where some basic facts concerning
curves in R d have been collected from [10].
relax the requirement that f should be differentiable but instead we constrain the length of
f . Note that by the definition of curves, f is still continuous. We give the following new
definition of principal curves.
curve f   is called a principal curve of length L for X if f   minimizes \Delta(f )
over all curves of length less than or equal to L.
A useful advantage of the new definition is that principal curves of length L always exist
if X has finite second moments, as the next result shows.
Assume that EkXk 2 ! 1. Then for any L ? 0 there exists a curve f   with
\Delta(f
The proof of the lemma is given in Appendix A.
Note that we have dropped the requirement of the HS definition that principal curves
be non-intersecting. In fact, Lemma 1 would not hold for non-intersecting curves of length
L without further restricting the distribution of X, since there are distributions for which
the minimum of \Delta(f ) is achieved only by an intersecting curve even though non-intersecting
curves can arbitrarily approach this minimum.
Remark: Connection with vector quantization Our new definition of principal curves
has been inspired by the notion of an optimal vector quantizer. The points y
are called the codepoints of an optimal k-point vector quantizer if
\Theta min
\Theta min

for any other collection of k points y In other words, the points y
k give
the best k-point representation of X in the mean squared sense. Optimal vector quantizers
are of great interest in lossy data compression, speech and image coding [11], and clustering
[12]. There is a strong connection between the definition of optimal vector quantizers and our
definition of a principal curve. Both minimize the same expected squared distance criterion,
while the vector quantizer is constrained to have at most k points, and we constrain the
length of a principal curve. This connection is further illuminated by a recent work of
Tarpey et al. [13] who define k points y to be self consistent if
are the Voronoi regions associated with y
kg. Thus our principal curves correspond to optimal vector quantizers
("principal points" by the terminology of [13]) while the HS principal curves correspond to
self consistent points.
While principal curves of a given length always exist, it appears difficult to demonstrate
concrete examples, unless the distribution of X is discrete or it is concentrated on a curve.
The same problem occurs in the theory of optimal vector quantizers, where except for the
scalar case the structure of optimal quantizers is unknown for even the most common
multivariate densities (e.g., see [11]).
Suppose now that n independent copies of X are given. These are called the
training data and they are assumed to be independent of X. The goal is to use the training
data to construct a curve of length at most L whose expected squared loss is close to that
of a principal curve for X.
Our method is based on a common model in statistical learning theory (e.g., see [9]).
We consider classes of curves increasing complexity. Given n data points
drawn independently from the distribution of X, we choose a curve as the estimator of the
principal curve from the kth model class S k by minimizing the empirical error. By choosing
the complexity of the model class appropriately as the size of the training data grows, the
chosen curve represents the principal curve with increasing accuracy.
We assume that the distribution of X is concentrated on a closed and bounded convex
set K ae R d . A basic property of convex sets in R d shows that there exists a principal curve
of length L inside K (see Lemma 2 in Appendix A), and so we will only consider curves in
K.
Let S denote the family of curves taking values in K and having length not greater than
L. For k - 1 let S k be the set of polygonal curves (broken lines) in K which have k segments
and whose lengths do not exceed L. Note that S k ae S for all k. Let
denote the squared distance between a point x 2 R d and the curve f . For any f 2 S the
empirical squared error of f on the training data is the sample average
where we have suppressed in the notation the dependence of \Delta n (f) on the training data. Let
our theoretical algorithm choose an f k;n 2 S k which minimizes the empirical error, i.e,
We measure the efficiency of f k;n in estimating f   by the difference J(f k;n ) between the
expected squared loss of f k;n and the optimal expected squared loss achieved by f   , i.e., we
let
\Delta(f
Our main result in this section proves that as the number
of data points n tends to infinity, and k is chosen to be proportional to n 1=3 , then J(f k;n )
tends to zero at a rate J(f k;n
Theorem 1 Assume that PfX 2 bounded and closed convex set K, let n
be the number of training points, and let k be chosen to be proportional to n 1=3 . Then the
expected squared loss of the empirically optimal polygonal line with k segments and length at
most L converges, as n !1, to the squared loss of the principal curve of length L at a rate
The proof of the theorem is given in Appendix B. To establish the result we use techniques
from statistical learning theory (e.g., see [14]). First, the approximating capability of the
class of curves S k is considered, and then the estimation (generalization) error is bounded via
covering the class of curves S k with ffl accuracy (in the squared distance sense) by a discrete
set of curves. When these two bounds are combined one obtains
r
where the term C(L; D; d) depends only on the dimension d, the length L, and the diameter
D of the support of X, but is independent of k and n. The two error terms are balanced by
choosing k to be proportional to n 1=3 which gives the convergence rate of Theorem 1.
Note that although the constant hidden in the O notation depends on the dimension d,
the exponent of n is dimension-free. This is not surprising in view of the fact that the class of
curves S is equivalent in a certain sense to the class of Lipschitz functions f
that Appendix A). It is known that the ffl-entropy, defined by the
logarithm of the ffl covering number, is roughly proportional to 1=ffl for such function classes
[15]. Using this result, the convergence rate O(n \Gamma1=3 ) can be obtained by considering ffl-covers
of S directly (without using the model classes S k ) and picking the empirically optimal curve
in this cover. However, the use of the classes S k has the advantage that they are directly
related to the practical implementation of the algorithm given in the next section.
3 A Polygonal Line Algorithm
Given a set of data points X the task of finding a polygonal line with
k segments and length L which minimizes 1
computationally difficult. We
propose a suboptimal method with reasonable complexity. The basic idea is to start with
a straight line segment f 1;n in each iteration of the algorithm we increase the
number of segments k by adding a new vertex to the polygonal line f k;n produced by the
previous iteration. After adding a new vertex, the positions of all vertices are updated in an
inner loop.
(a)
-0.4
(b)
-0.4
(c)
-0.4
(d)
-0.4

Figure

2: The curves f k;n produced by the polygonal line algorithm for data points. The
data was generated by adding independent Gaussian errors to both coordinates of a point chosen
randomly on a half circle. (a) f 1;n , (b) f 2;n , (c) f 4;n , (d) f 11;n (the output of the algorithm).
The inner loop consists of a projection step and an optimization step. In the projection
step the data points are partitioned into "nearest neighbor regions" according to which
segment or vertex they project. In the optimization step the new position of a vertex v i is
determined by minimizing an average squared distance criterion penalized by a measure of
the local curvature, while all other vertices are kept fixed. These two steps are iterated so
that the optimization step is applied to each vertex v i , in a cyclic fashion
(so that after v k+1 , the procedure starts again with v 1 ), until convergence is achieved and
f k;n is produced. Then a new vertex is added.
The algorithm stops when k exceeds a threshold c(n; \Delta). This stopping criterion is based
on a heuristic model complexity measure, determined by the number segments k, the number
of data points n, and the average squared distance \Delta n (f k;n ).
The flow-chart of the algorithm is given in Figure 3. The evolution of the curve produced
by the algorithm is illustrated in Figure 2. We note here that the objective function to be
minimized in the vertex optimization procedure is based partly on heuristic considerations.
As explained in Section 3.3, the algorithm in this step searches for a (local) minimum of a
the average squared distance penalized by the local curvature. The heuristic lies in the data
dependent form of the penalty factor. Similarly to the HS algorithm, we have no formal
proof that the practical algorithm will converge, but in practice, after extensive testing, it is
observed to converge.
Convergence?
Y
Y
END
Projection
Initialization
Add new vertex
Vertex optimization

Figure

3: The flow chart of the polygonal line algorithm.
3.1 The Initialization Step
To obtain f 1;n , we take the shortest segment of the first principal component line which
contains all of the projected data points. To keep the computational complexity low, we
compute the first principal component of a constant number of points randomly chosen from
the n data points. This choice suffices for our purposes since the algorithm needs only a
reasonable approximation of the first principal component.
3.2 The Projection Step
Let f denote a polygonal line with vertices closed line segments s
such that s i connects vertices v i and v i+1 . In this step the data set X n is partitioned into
(at the nearest neighbor regions of
the vertices and segments of f , in the following manner. For any x 2 R d let \Delta(x; s i ) be the
squared distance from x to s i (see definition (3)), and let \Delta(x; . Then we let
Upon setting
are defined by
The resulting partition is illustrated in Figure 4.001101
s
s
s

Figure

4: The nearest neighbor partition of R 2 induced by the vertices and segments of f .
3.3 The Vertex Optimization Step
In this step the new position of a vertex v i is determined. In the theoretical algorithm
the average squared distance \Delta n (x; f) is minimized subject to the constraint that f is a
polygonal line with k segments and length not exceeding L. One could use a Lagrangian
formulation and attempt to find a new position for v i (while all other vertices are fixed) such
that the penalized squared error \Delta is minimum. However, we have observed
that this approach is very sensitive to the choice of -. On the other hand, most principal
curve applications require a smooth curve solution. To avoid over-fitting, HS used scatterplot
or spline smoothing. We chose to penalize the local curvature to obtain smoother curves.
Due to the fact that only one vertex is moved at a time, penalizing the curvature will also
implicitly penalize the length of the curve. After considering several possibilities, we found
that the following measures of local curvature work especially well. At inner vertices v i
we penalize the sum of the cosines of the three angles at vertices v
At the endpoints and at their immediate neighbors (v i , 1), the
penalty on a nonexistent angle is replaced by the squared length of the first (or last) segment.
Formally, let fl i denote the angle at vertex v i , let -(v
and let . Then the penalty P (v i ) at vertex v i is given by
The local measure of the average squared distance is calculated from the data points which
project to v i or to the line segment(s) starting at v i (see Projection Step). Accordingly, let
define the local average squared distance as a function of v i by
We use a gradient (steepest descent) method to minimize
This part of the algorithm is modular, i.e., the simple procedure we are
using can be substituted with a more sophisticated nonlinear programming procedure at the
expense of increased computational complexity.
One important issue is the amount of smoothing required for a given data set. In the HS
algorithm one needs to determine the penalty coefficient of the spline smoother, or the span
of the scatterplot smoother. In our algorithm, the corresponding parameter is the curvature
penalty factor - p . If some a priori knowledge about the distribution is available, one can
use it to determine the smoothing parameter. However in the absence of such knowledge,
the coefficient should be data-dependent. Intuitively, - p should increase with the number
of segments and the size of the average squared error, and it should decrease with the data
size. Based on heuristic considerations and after carrying out practical experiments, we set
p is a parameter of the algorithm.
3.4 Adding a New Vertex
We start with the optimized f k;n and choose the segment that has the largest number of
data points projecting to it. The midpoint of this segment is selected as the new vertex.
Formally, let j. Then the new vertex is v new
Stopping Condition
According to the theoretical results of Section 2, the number of segments k should be proportional
to n 1=3 to achieve the O(n 1=3 ) convergence rate for the expected squared dis-
tance. Though the theoretical bounds are not tight enough to determine the optimal number
of segments for a given data size, we found that k - n 1=3 also works in practice.
To achieve robustness we need to make k sensitive to the average squared distance. The
stopping condition blends these two considerations. The algorithm stops when k exceeds
3.6 Computational Complexity
The complexity of the inner loop is dominated by the complexity of the projection step, which
is O(nk). Increasing the number of segments by one at a time (as described in Section 3.4),
the complexity of the algorithm to obtain f k;n is O(nk 2 ). Using the stopping condition of
Section 3.5, the computational complexity of the algorithm becomes O(n 5=6 ). This is slightly
better than the O(n 2 ) complexity of the HS algorithm.
The complexity can be dramatically decreased in certain situations. One possibility is to
add more than one vertex at a time. For example, if instead of adding only one vertex, a new
vertex is placed at the midpoint of every segment, then we can reduce the computational
complexity for producing f k;n to O(nk log k). One can also set k to be a constant if the
data size is large, since increasing k beyond a certain threshold brings only diminishing
returns. These simplifications work well in certain situations, but the original algorithm is
more robust.
4 Experimental Results
We have extensively tested our algorithm on two-dimensional data sets. In most experiments
the data was generated by a commonly used (see, e.g., [2] [6] [7]) additive model
where Y is uniformly distributed on a smooth planar curve (hereafter called the generating
curve) and e is bivariate additive noise which is independent of Y.
Since the "true" principal curve is not known (note that the generating curve in the model
e is in general not a principal curve either in the HS sense or in our definition),
it is hard to give an objective measure of performance. For this reason, in what follows, the
performance is judged subjectively, mainly on the basis of how closely the resulting curve
follows the shape of the generating curve.
In general, in simulation examples considered by HS the performance of the new algorithm
is comparable with the HS algorithm. Due to the data-dependence of the curvature penalty
factor and the stopping condition, our algorithm turns out to be more robust to alterations
in the data generating model, as well as to changes in the parameters of the particular model.
We use varying generating shapes, noise parameters, and data sizes to demonstrate the
robustness of the polygonal line algorithm. All plots show the generating curve (Generator
Curve), the curve produced by our polygonal line algorithm (Principal Curve), and the curve
produced by the HS algorithm with spline smoothing (HS Principal Curve), which we have
found to perform better than the HS algorithm using scatterplot smoothing. For closed
generating curves we also include the curve produced by the Banfield and Raftery (BR)
algorithm [5], which extends the HS algorithm to closed curves (BR Principal Curve). The
two coefficients of the polygonal line algorithm are set in all experiments to the constant
values
plots have been normalized to fit in a 2 \Theta 2 square. The
parameters given below refer to values before this normalization.
In

Figure

5 the generating curve is a circle of radius
bivariate uncorrelated Gaussian with variance E(e 2
2. The performance of
the three algorithms (HS, BR, and the polygonal line algorithm) is comparable, although
the HS algorithm exhibits more bias than the other two. Note that the BR algorithm [5] has
been tailored to fit closed curves and to reduce the model bias. In Figure 6, only half of the
circle is used as a generating curve and the other parameters remain the same. Here, too,
both the HS and our algorithm behave similarly.
When we depart from these usual settings the polygonal line algorithm exhibits better
behavior than the HS algorithm. In Figure 7(a) the data set of Figure 6 was linearly transformed
using the matrix ( 0:6 0:6
\Gamma1:0 1:2 ). In Figure 7(b) the transformation
1:0 \Gamma0:2
was used.
The original data set was generated by an S-shaped generating curve, consisting of two half
circles of unit radii, to which the same Gaussian noise was added as in Figure 6. In both
cases the polygonal line algorithm produces curves that fit the generator curve more closely.
This is especially noticeable in Figure 7(a) where the HS principal curve fails to follow the
shape of the distorted half circle.
There are two situations when we expect our algorithm to perform particularly well. If the
distribution is concentrated on a curve, then according to both the HS and our definitions the
principal curve is the generating curve itself. Thus, if the noise variance is small, we expect
both algorithms to very closely approximate the generating curve. The data in Figure 8(a)
was generated using the same additive Gaussian model as in Figure 5, but the noise variance
was reduced to E(e 2
2. In this case we found that the polygonal line
algorithm outperformed both the HS and the BR algorithms
The second case is when the sample size is large. Although the generating curve is not
necessarily the principal curve of the distribution, it is natural to expect the algorithm to
well approximate the generating curve as the sample size grows. Such a case is shown in

Figure

8(b), where data points were generated (but only a small subset of these
was actually plotted). Here the polygonal line algorithm approximates the generating curve
with much better accuracy than the HS algorithm.
5 Conclusion
A new definition of principal curves has been offered. The new definition has significant
theoretical appeal; the existence of principal curves under this definition can be proved
under very general conditions, and a learning method for constructing principal curves for
finite data sets yields to theoretically analysis.
Inspired by the new definition and the theoretical learning scheme, we have introduced
a new practical polygonal line algorithm for designing principal curves. Lacking theoretical
results concerning both the HS and our polygonal line algorithm, we compared the two
methods through simulations. We have found that in general our algorithm has performance
either comparable with the performance of the original HS algorithm and it exhibits better,
more robust behavior when the data generating model is varied. It should be mentioned that
these findings cannot be called entirely conclusive due mainly to the absence of an objective
performance measure. In practical applications, each method can have different advantages.
In this respect, we believe that the new principal curve algorithm may prove useful where a
compact and accurate description of a pattern or an image is required, e.g., in skeletonization
of handwritten characters or in feature extraction. These are issues for future work.


Appendix

A
Curves in R d
continuous mapping (curve). The length of f over an interval
denoted by l(f ; ff; fi), is defined by
where the supremum is taken over all finite partitions of [ff; fi] with arbitrary subdivision
points 1. The length of f over its entire domain
[a; b] is denoted by l(f ). If l(f) ! 1, then f is said to be rectifiable. It is well known
that rectifiable iff each coordinate function R is of bounded
variation.
Two curves are said to be equivalent if there exist two
nondecreasing continuous real functions OE
In this case we write f - g, and it is easy to see that - is an equivalence relation. If
curve g over [a; b] is said to be parametrized by its arc length if
a for any a - t - b. Let f be a curve over [a; b] with length L. It is not hard
to see that there exists a unique arc length parametrized curve g over [0; L] such that f - g.
Let f be any curve with length L 0 - L, and consider the arc length parametrized curve
- f with parameter interval [0; L 0 ]. By definition (A.1), for all s
g, and -
satisfies
the following Lipschitz condition: For all
On the other hand, note that if -
g is a curve over [0; 1] which satisfies the Lipschitz condition
(A.2), then its length is at most L.
Let f be a curve over [a; b] and denote the squared Euclidean distance from any x 2 R d
to f by
a-t-b
Note that if l(f) ! 1, then by the continuity of f , its graph
is a compact subset of R d , and the infimum above is achieved for some t. Also, since G f
if f - g, we also have that \Delta(x; f) = \Delta(x; g) for all g - f .
Proof of Lemma 1 Define
First we show that the above infimum does not change if we add the restriction that all f
lie inside a closed sphere of large enough radius r and centered at
the origin. Indeed, without excluding nontrivial cases, we can assume that \Delta   ! EkXk 2 .
Denote the distribution of X by - and choose r ? 3L large enough such that
Z
for some ffl ? 0. If f is such that G f
is not entirely contained in S(r), then for all x 2 S(r=3)
we have \Delta(x; f) ? kxk 2 since the diameter of G f is at most L. Then (A.3) implies that
\Delta(f
Z
and thus
ae S(r)g: (A.4)
In view of (A.4) there exists a sequence of curves ff n g such that l(f n
ae S(r)
for all n, and \Delta(f n . By the discussion preceding (A.2), we can assume without loss
of generality that all f n are defined over [0; 1] and
Consider the set of all curves C over [0; 1] such that f 2 C iff kf(t 1
ae S(r). It is easy to see that C is a closed
set under the uniform metric d(f ; Also, C is an equicontinuous
family of functions and sup t kf(t)k is uniformly bounded over C. Thus C is a compact metric
space by the Arzela-Ascoli theorem (see, e.g., [16]). Since f n 2 C for all n, it follows that
there exists a subsequence f n k
converging uniformly to an f   2 C.
To simplify the notation let us rename ff n k
g as ff n g. Fix x 2 R d , assume \Delta(x; f n
\Delta(x; f   ), and let t x be such that \Delta(x; f   Then by the triangle inequality,
)k:
By symmetry, a similar inequality holds if \Delta(x; f n
EkXk 2 is finite, there exists A ? 0 such that
and therefore
Since the Lipschitz condition on f   guarantees that l(f   ) - L, the proof is complete.
Assume that PfX 2 closed and convex set K, and let f be a curve
with l(f) - L. Then there exists a curve - f such that G - f
ae K, l( - f) - L, and
Proof. For each t in the domain of f , let - f(t) be the unique point in K such that
It is well known that - f(t) satisfies
where h\Delta; \Deltai denotes the usual inner product in R d (see, e.g, [17]). Then for all
where the inequality follows from (A.6) since - f(t 1 continuous (it is
a curve) and similar inequality shows that for all t and x 2 K,
so that \Delta( - f


Appendix

Proof of Theorem 1. Let f
k denote the curve in S k minimizing the squared loss, i.e.,
f
\Delta(f
The existence of a minimizing f
k can easily be shown using a simpler version of the proof of
Lemma 1. Then J(f k;n ) can be decomposed as
where, using standard terminology, \Delta(f k;n
k ) is called the estimation error and \Delta(f
\Delta(f   ) is called the approximation error. We consider these terms separately first, and then
choose k as a function of the training data size n to balance the obtained upper bounds in
an asymptotically optimal way.
Approximation
For any two curves f and g of finite length define their (nonsymmetric) distance by
min
s
Note that ae( - f ; -
- g, i.e., ae(f ; g) is independent of the particular
choice of the parametrization within equivalence classes. Next we observe that if the diameter
of K is D, and G f ; G g 2 K, then for all x 2 K,
and therefore
To prove (B.1), let x 2 K and choose t 0 and s 0 such that \Delta(x; f)
be an arbitrary arc length parametrized curve over [0; L 0 ], where L 0 - L.
as a polygonal curve with vertices
some s, we have
min
s
Note that l(g) - L 0 , by construction, and thus g 2 S k . Thus for every f 2 S there exists a
such that ae(f ; g) - L=(2k). Now let g 2 S k be such that ae(f   ; g) - L=(2k). Then by
(B.2) we conclude that the approximation error is upper bounded as
\Delta(f
Estimation
For each ffl ? 0 and k - 1 let S k;ffl be a finite set of curves in K which form an ffl-cover of
S k in the following sense. For any f 2 S k there is an f 0 2 S k;ffl which satisfies
sup
The explicit construction of S k;ffl is given in Appendix C. Since f k;n 2 S k (see (5)), there
exists an f 0
We introduce the
compact notation X for the training data. Thus we can write
where (B.5) follows from the approximating property of f 0
k;n and the fact that the distribution
of X is concentrated on K. (B.6) holds because f k;n minimizes \Delta n (f) over all f 2 S k , and
is an ordinary expectation
of the type E [\Delta(X; f )], f 2 S k;ffl . Thus for any t ? 2ffl the union bound implies
where jS k;ffl j denotes the cardinality of S k;ffl .
Recall now Hoeffding's inequality [18] which states that if Y are independent
and identically distributed real random variables such that 0 - Y i - A with probability one,
then for all u ? 0,
Since the diameter of K is D, we have such that
K. Thus 0 - \Delta(X; f) - D 2 with probability one and by Hoeffding's inequality for all
we have
which implies by (B.8) that
for any t ? 2ffl. Using the identity
nonnegative random
variable Y , we can write for any u ? 0,
\Delta(f k;n
Z 1P
dt
r
where (B.10) follows from the inequality
follows by setting
, where log denotes natural logarithm. The following
lemma, which is proved in Appendix C, demonstrates the existence of a suitable covering
set S k;ffl .
Lemma 3 For any ffl ? 0 there exists a finite collection of curves S k;ffl in K such that
sup
and
d
d
d
LD
d
d
where V d is the volume of the d-dimensional unit sphere and D is the diameter of K.
It is not hard to see that setting gives the upper bound
log(jS
where C(L; D; d) does not depend on k. Combining this with (B.11) and the approximation
bound given by (B.3) results in
\Delta(f k;n
r
The rate at which \Delta(f k;n ) approaches \Delta(f   ) is optimized by setting the number of segments
k to be proportional to n 1=3 . With this choice J(f k;n has the asymptotic
convergence rate
and the proof of Theorem 1 is complete.


Appendix

Proof of Lemma 3 Consider a rectangular grid with side length ffi ? 0 in R d . With each
point y of this grid associate its Voronoi region (a hypercube of side length ffi), defined as the
set of points which are closer to y than to any other points of the grid. Let K ffi ae K denote
the collection of points of this grid which fall in K plus the projections of those points of the
grid to K whose Voronoi regions have a nonempty intersections with K. Then we clearly
have
min
dffi: (C.1)
d) and define S k;ffl to be the family of all polygonal curves - f having k
vertices
satisfying the length constraint
dffi: (C.2)
To see that S k;ffl has the desired covering property, let f 2 S k be arbitrary with vertices
f be the
polygonal curve with vertices -
by the definition of S k ,
the triangle inequality implies that - f satisfies (C.2) and thus - f 2 S k;ffl . On the other hand,
without loss of generality assume that the line segment connecting y i\Gamma1 and y i and the line
segment connecting -
y
y i are both linearly parametrized over [0; 1]. Then
dffi:
This shows that
dffi=2. Then it follows from (B.1) that S k;ffl is an
ffl-cover for S k since for all x 2 K,
denote the length of the ith segment of - f and let
where dxe denotes the least integer not less than x. Fix the sequence -
define S
as the set of all - f 2 S k;ffl whose segment lengths generate this particular
sequence. To bound jS k;ffl ( -
note that the first vertex -
y 0 of an - f 2 S k;ffl ( -
L can be any of
the points in K ffi which contains as many points as there are Voronoi cells intersecting K.
Since the diameter of K is D, there exists a sphere of radius D
which contains these
Voronoi cells. Thus the cardinality of K ffi can be upper bounded as
where V d is the volume of the unit sphere in R d . Assume -
y has been
chosen. Since k- y
there are no more than
possibilities for choosing -
y i . Therefore,
d
Y
By (C.2) and the definition of the -
Therefore the arithmetic-geometric mean inequality implies that
d
Y
and thus
d
d
On the other hand, by (C.3) we have
+2k and therefore the number of distinct
sequences -
L kis upper bounded by
Substituting
d) we obtain
d
d
d
LD
d
d



--R

"Principal curves and surfaces."
"Principal curves,"
"An algorithm for vector quantizer design,"
"Robust locally weighted regression and smoothing scatterplots,"
"Ice floe identification in satellite images using mathematical morphology and clustering about principal curves,"
"Principal curves revisited,"
"Self-organization as an iterative kernel smoothing pro- cess,"
"Another look at principal curves and surfaces."
The Nature of Statistical Learning Theory.
Introductory Real Analysis.

Clustering Algorithms.
"Principal points and self-consistent points of elliptical distributions,"
A Probabilistic Theory of Pattern Recognition.
"ffl-entropy and ffl-capacity of sets in function spaces,"
Real analysis and probability.
Optimization by vector space methods.
"Probability inequalities for sums of bounded random variables,"
--TR

--CTR
Balzs Kgl , Adam Krzyak, Piecewise Linear Skeletonization Using Principal Curves, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.1, p.59-74, January 2002
Peter Meinicke , Stefan Klanke , Roland Memisevic , Helge Ritter, Principal Surfaces from Unsupervised Kernel Regression, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.9, p.1379-1391, September 2005
B. Bhushan , J. A. Romagnoli, A strategy for feature extraction of high dimensional noisy data, Proceedings of the 25th IASTED international conference on Modeling, indentification, and control, p.441-445, February 06-08, 2006, Lanzarote, Spain
Zhiguo Cheng , Mang Chen , Yuncai Liu, A robust algorithm for image principal curve detection, Pattern Recognition Letters, v.25 n.11, p.1303-1313, August 2004
J. J. Verbeek , N. Vlassis , B. Krse, A k-segments algorithm for finding principal curves, Pattern Recognition Letters, v.23 n.8, p.1009-1017, June 2002
Jos Koetsier , Ying Han , Colin Fyfe, Twinned principal curves, Neural Networks, v.17 n.3, p.399-409, April 2004
B. S. Y. Lam , H. Yan, A curve tracing algorithm using level set based affine transform, Pattern Recognition Letters, v.28 n.2, p.181-196, January, 2007
Hujun Yin, Data visualisation and manifold mapping using the ViSOM, Neural Networks, v.15 n.8-9, p.1005-1016, October 2002
Jochen Einbeck , Gerhard Tutz , Ludger Evers, Local principal curves, Statistics and Computing, v.15 n.4, p.301-313, October   2005
Alexander J. Smola , Sebastian Mika , Bernhard Schlkopf , Robert C. Williamson, Regularized principal manifolds, The Journal of Machine Learning Research, 1, p.179-209, 9/1/2001
Kui-Yu Chang , J. Ghosh, A Unified Model for Probabilistic Principal Surfaces, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.1, p.22-41, January 2001
curve fitting;feature extraction;learning systems;unsupervised learning;piecewise linear approximation;vector quantization