Inference on covariance operators via concentration inequalities:
k-sample tests, classification, and clustering via Rademacher
complexities

arXiv:1604.06310v1 [stat.ME] 21 Apr 2016

Adam B. Kashlak1 , John A. D. Aston 2 , and Richard Nickl2
1

Cambridge Centre for Analysis, University of Cambridge, Wilberforce Rd, Cambridge, CB3 0WB, U.K
2
Statistical Laboratory, University of Cambridge, Wilberforce Rd, Cambridge, CB3 0WB, U.K

March 27, 2018
Abstract
We propose a novel approach to the analysis of covariance operators making use of concentration inequalities. First,
non-asymptotic confidence sets are constructed for such operators. Then, subsequent applications including a k sample
test for equality of covariance, a functional data classifier, and an expectation-maximization style clustering algorithm are
derived and tested on both simulated and phoneme data.

1

Introduction

Functional data spans many realms of applications from medical imaging, (Jiang et al., 2014), to speech and linguistics,
(Pigoli et al., 2014), to genetics, (Panaretos et al., 2010). General inference techniques for functional data are one area of
analysis that has received much attention in recent years from the construction of confidence sets, to other topics such as
k-sample tests, classification, and clustering of functional data. Most testing methodology treats the data as continuous L2
valued functions and subsequently reduces the problem to a finite dimensional one through expansion in some orthogonal
basis such as the often utilized Karhunen-Loève expansion (Horváth and Kokoszka, 2012). However, inference in the more
general Banach space setting, making use of a variety of norms, has not been addressed adequately. We propose a novel
methodology for performing fully functional inference through the application of concentration inequalities; for general
concentration of measure results, see Ledoux (2001) and Boucheron et al. (2013). Special emphasis is given to inference on
covariance operators, which offers a fruitful way to analyze functional data.
Imagine multiple samples of speech data collected from multiple speakers. Each speaker will have his or her own
sample covariance operator taking into account the unique variations of his or her speech and language. An exploratory
researcher may want to find natural clusters amidst the speakers perhaps corresponding to gender, language, or regional
accent. Meanwhile, a linguist studying the similarities between languages may want to test for the equality of such
covariances. A computer scientist may need to implement an algorithm that when given speech data quickly identifies what
language is being spoken and furthermore parses the sound clip and identifies each individual phoneme in order to process
the speech into text. Our proposed method has the versatility to yield statistical tests that answer these questions as well as
others.
Past methods for analyzing covariance operators (Panaretos et al., 2010; Fremdt et al., 2013) rely on the Hilbert-Schmidt
setting for their inference. However, the recent work of Pigoli et al. (2014) argues that the use of the Hilbert-Schmidt
metric ignores the geometry of the covariance operators and that more statistical power can be gained by using alternative
metrics. Hence, we approach such inference for covariance operators in the full generality of Banach spaces, by using
a non-asymptotic concentration of measure approach, which has previously been used in nonparametric statistics and
machine learning, sometimes under the name of ‘Rademacher complexities’ (Koltchinskii, 2001, 2006; Bartlett et al., 2002;
Bartlett and Mendelson, 2003; Giné and Nickl, 2010; Arlot et al., 2010; Lounici and Nickl, 2011; Kerkyacharian et al.,
2012; Fan, 2011). These concentration inequalities provide a natural way to construct non-asymptotic confidence regions.

1

2

Definitions and notation

Generally, we will consider functional data to be in the Hilbert space L2 (I) for I ⊂ R. Our estimated covariance operators
of interest are treated as Banach space valued random variables. Let

Op(L2 ) = T : L2 (I) → L2 (I) ∃M ≥ 0, kT φkL2 ≤ M kφkL2 ∀φ ∈ L2 (I)
denote the space of all bounded linear operators mapping L2 into L2 . This is where our covariance operators will live.
The metrics that will be investigated are those that correspond to the p-Schatten norms. When p 6= 2, these are not
Hilbert norms. Hence, we will consider the general Banach space setting.
Definition 2.1 (p-Schatten Norm) Given separable Hilbert spaces H1 and H2 , a bounded linear operator Σ : H1 → H2 ,
p
and p ∈ [1, ∞), then the p-Schatten norm is kΣkp = tr((Σ∗ Σ)p/2 ). For p = ∞, the Schatten norm is the operator norm:
kΣk∞ = supf ∈H1 (kΣf kH2 /kf kH1 ). In the case that Σ is compact, self-adjoint, and trace-class, then given the associated
p
eigenvalues {λi }∞
i=1 , the p-Schatten norm coincides with the ` norm of the eigenvalues:

P∞
p
p
kλk`p = i=1 |λi | , p ∈ [1, ∞)
p
kΣkp =
maxi∈N |λi | ,
p=∞
In order to construct a covariance operator from a sample of functional data, the notion of tensor product is required.
Let f, g ∈ L2 (I) and φ in the dual space L2 (I)∗ with inner product hf, φi = φ(f ). The tensor product, f ⊗ g, is the rank
one operator defined by (f ⊗ g)φ = hg, φi f = φ(g)f .
Secondly, we will implement a Rademacher symmetrization technique in the concentration inequalities to come. This
requires the use of the namesake Rademacher random variables.
Definition 2.2 (Rademacher Distribution) A random variable ε ∈ R has a Rademacher distribution if pr (ε = 1) =
pr (ε = −1) = 1/2.
One particularly fruitful avenue of functional data analysis is the analysis of covariance operators. Such an approach to
functional data has been discussed by Panaretos et al. (2010) for DNA microcircles, by Fremdt et al. (2013) for the egg
laying practices of fruit flies, and by Pigoli et al. (2014) with application to differentiating spoken languages.


2
Definition 2.3 (Covariance Operator) Let I ⊆ R, and let f be a random function (variable) in L2 (I) with E kf kL2 <

∞ and mean zero. The associated covariance operator Σf ∈ Op(L2 ) is defined as Σf = E f ⊗2 = E (hf, ·i f ) .
m
If I = {i1 , . . . , im } has finite cardinality, then f = (f
 1 , . . . , fm ) is a random vector in R and for some fixed vector
v ∈ Rm , E (hf, vi f ) = E f f T v where Σf = E f f T is the usual covariance matrix. Covariance operators are integral
operators with the kernel function cf (s, t) = cov{f (s), f (t)} ∈ L2 (I × I). Such operators are characterized by the result
that for f ∈ L2 (I), Σf is a covariance operator if and only if it is trace-class, self-adjoint, and compact on L2 (I) where the
symmetry follows immediately from the definition and the finite trace norm comes from Parseval’s equality (Bosq, 2012,
Theorem 1.7)(Horváth and Kokoszka, 2012, Section 2.3).

4
Furthermore, working under the assumption that E kf kL2 < ∞, we will require tensor powers of covariance
2
operators denoted as Σ⊗2 : Op(L2 ) → Op(L2 ). For a basis {ei }∞
⊗ ej }∞
i=1 ∈ L (I) with corresponding basis {eiP
i,j=1 for
∞
2
⊗2
2
Op(L (I)), the previous definition is extended to Σ = hΣ, ·i Σ where for Σ1 , Σ2 ∈ Op(L ) with Σ1 = i,j=1 λi,j ei,j
P∞
P∞
and Σ2 = i,j=1 γi,j ei,j , then hΣ1 , Σ2 i = i,j λi,j γi,j . Specifically for covariance operators, the tensor power takes on
a similar integral operator form with kernel cΣ (s, t, u, v) = cov (f (s), f (t)) cov (f (u), f (v)) .
Given an Hilbert space H with inner product h·, ·i, the adjoint of a bounded linear operator Σ : H → H, denoted
as Σ∗ , is the unique operator such that hΣf, gi = hf, Σ∗ gi for f, g ∈ H. The existence of which is given by the Riesz
representation theorem (Rudin, 1987, chapter 2). For self-adjoint operators, such as the covariance operators of interest,
Σ = Σ∗ .
We begin with a sample of functional data. Let f1 , . . . , fn ∈ L2 (I) be independent and
Pn identically distributed
observations with mean zero and covariance operator Σ. Let the sample mean be f¯ = n−1 i=1 fi and the empirical
Pn
estimate of Σ be Σ̂ = n−1 i=1 (fi − f¯) ⊗ (fi − f¯). The initial goal is to construct a confidence set for Σ with respect to
some metric d(·, ·) of the form {Σ : d(Σ̂, Σ) ≤ r(n, Σ̂, α)}, which has coverage 1 − α for any desired α ∈ [0, 1] and a
radius r depending only on the data and α. Such a confidence set can be utilized for a wide variety of statistical analyses.

2

3
3.1

Constructing a confidence set
Confidence sets for the mean in Banach spaces

The goal of this section is to construct a non-asymptotic confidence region in the Banach space setting. These can
subsequently be specialized to our case of interest, covariance operators, when the Xi below are replaced with fi⊗2 .
The construction of our confidence set is based on Talagrand’s concentration inequality (Talagrand, 1996) with
explicit constants. This inequality is typically stated for empirical processes (Giné and Nickl, 2016, Theorem 3.3.9 and
3.3.10 ), but applies to random variables with values in a separable Banach space (B, k·kB ) as well by simple duality
arguments (Giné and Nickl, 2016, Example 2.1.6). Let X1 , . . . , Xn ∈ (B, k·kB ) be mean zero independent and identically
distributed Banach space valued random variables with kXi kB ≤ U for all i = 1, . . . , n where U is some positive constant.
Furthermore, let h·, ·i : B × B ∗ → R such that for X ∈ B and φ ∈ B ∗ then hX, φi = φ(X). Define
Z=

sup

n
X

n
X

hXi , φi =

kφkB ∗ ≤1 i=1

n

Xi

i=1

σ2 =

,
B



1X
2
sup E hXi , φi ,
n i=1 kφkB∗ ≤1

where the supremum is taken over a countably dense subset of the unit ball of B ∗ . Furthermore,
define v = 2U E (Z) + nσ 2 .

2
Then, pr{Z > E (Z) + r} ≤ exp{−r /(2v + 2rU/3)}. Rewriting Z as n X̄ − E X̄ B results in







−n2 r2
pr X̄ − E X̄ B > E X̄ − E X̄ B + r < exp
2v + 2nrU/3

where kXi kB < U and v = 2nU E{ X̄ − E X̄ B } + nσ 2 .

The above tail bound incorporates the unknown E(kX̄ − E X̄ kB ). Consequently,
a symmetrization technique is used.
Pn
This term is replaced by the norm of the Rademacher average Rn = n−1 i=1 εi (Xi − X̄) where the εi are independent
and identically distributed Rademacher random variables also independent of the Xi . This substitution is justified by
invoking the symmetrization inequality (Giné and Nickl, 2016, Theorem 3.1.21),
)
(
)
(
n
n

1X
1X
≤ 2E
= 2E (kRn kB ) .
E (Z) = E
(Xi − E X̄ )
εi (Xi − X̄)
n i=1
n i=1
B

B

In practice, the coefficient of 2 is unnecessary as averaging the data already has a symmetrizing effect, and if starting
from symmetric data, then X − X̄ and ε(X − X̄) are equal in distribution. This result allows us to replace the original
expectation with the expectation of the Rademacher average. Furthermore, Talagrand’s inequality also applies to Rn .
Hence, the Rademacher average concentrates strongly about its expectation, which justifies dropping the expectation. In
practice, one can use the intermediary Eε (kRn kB ), which can be approximated for reasonable sized data sets via Monte
Carlo simulations of the εi . However, this is not strictly necessary, and for large data sets, a single random draw of εi will
suffice (Giné and Nickl, 2016, Section 3.4.2).
The resulting (1 − α)-confidence set is
(
)


 1/2 U log(2α)
2
2
X : X − X̄ B ≤ kRn kB +
log(2α) σ + 2U kRn kB
+
.
(3.1)
n
3n
To make use of these results in practice, both the weak variance σ 2 must be estimated for the data and a reasonable choice
of U must be made, and a main contribution of this present paper is to propose some theoretically motivated but practically
useful non-asymptotic choices for these constants that work for the functional data applications we are investigating. This
will be discussed in the setting of covariance operators in the following subsection.

3.2

Confidence sets for covariance operators

To construct a confidence set for covariance operators, let our functional data fi ∈ L2 (I) and fi⊗2 = fi ⊗ fi ∈ Op(L2 ),
the Hilbert space of bounded linear operators mapping L2 to L2 , such that (fi ⊗ fi )φ = hf, φi f for some φ ∈ L2 . Hence,
for some desired p-Schatten norm, k·kp , with p ∈ [1, ∞) and with conjugate q = p/(p − 1), we have
n

Z=

1X
fi ⊗ fi − E (fi ⊗ fi )
n i=1

= Σ̂ − Σ

σ2 =

,
p

p

3

n
n
o

1X
2
sup E fi⊗2 − E fi⊗2 , Π
n i=1 kΠkq ≤1

for the supremum
over a countably dense subset of the unit ball of Π ∈ Op(L2 ). The Rademacher average,
Pn being taken
−1
⊗2
¯
Rn = n
− Σ̂}, is also in Op(L2 ), because for any φ ∈ L2 (I) and for some M ∈ R, kRn φkL2 ≤
i=1 εi {(fi − f )
|εi | k{(fi − f¯)⊗2 − Σ̂}φkL2 ≤ M kφkL2 since {(fi − f¯)⊗2 − Σ̂} is a bounded operator and |εi | = 1. Furthermore, in
the case that there exists a fixed c ∈ R with kfi kL2 ≤ c for all i corresponding to a physical bound on the energy of
fi , kfi⊗2 kp = kfi k2L2 ≤ c2 = U . It will be determined below that U ≥ σ in this case. In general, setting U = σ gives
good experimental results when fi is Gaussian as will be discussed in later sections. This results in vn ≈ σ 2 /n. For any
p ∈ [1, ∞) and α ∈ [0, 1/2], the proposed (1 − α)-confidence set inspired by Equation 3.1 for covariance operators is
"
#

1/2
2
σ log(2α)
Cn,1−α = Σ : Σ̂ − Σ ≤ kRn kp + σ − log(2α)
−
.
n
3n
p
where σ depends on the distribution on the functional data. As a rule of thumb for the choice of σ 2 , as we will now show, is
to note that σ 2 ≤ kE(f ⊗4 ) − Σ⊗2 kp and to estimate this bound empirically by σ̂ 2 . For example, when the fi are from a
Gaussian process σ̂ ≤ 21/2 kΣkp as explained in detail in Section 3.4.
To calculate the weak variance σ 2 , define f ⊗n = f ⊗ . . . ⊗ f to be the n-fold tensor product of f with itself and
2
4
extend the definition of h·, ·i : (L2 )⊗4 × {(L2 )⊗4 }∗ → R such that f ⊗4 , φ⊗4 = f ⊗2 , φ⊗2 = hf, φi = φ(f )4 . For
operators Π ∈ {(L2 )⊗2 }∗ and Ξ ∈ {(L2 )⊗4 }∗ , the weak variance is
σ2 =

n
n
o

1X
2
sup E fi⊗2 − E fi⊗2 , Π
n i=1 kΠkq ≤1
n

≤

D
 

1X
sup E fi⊗4 − E f ⊗2
n i=1 kΞkq ≤1

⊗2

E

, Ξ ≤ E f ⊗4 − Σ⊗2

p

where the inequality stems from the fact that the supremum is being taken over a larger set. However, in the Hilbert space
setting, the dual of the tensor product does coincide with the tensor product of the dual space, and thus the above inequality
2
can be replaced with an equality if the Hilbert-Schmidt norm, 2-Schatten norm, is used. Given a bound kfi kL2 ≤ c2 = U ,
2
⊗4
4
4
2
then σ ≤ kE f
kp ≤ E(kf kL2 ) ≤ c = U .

3.3

The weak variance for p = ∞

Let E be a countable dense subset of the unit ball of L2 (I). In the case p = ∞, we cannot use duality, but can still write Z
and σ 2 as suprema over the countable set and achieve the same results as above.
Z=
σ2 =
≤
As before, if fi⊗2

3.4

n
D
E
X
 ⊗2

1
fi − E fi⊗2 e, e = sup (Σ̂ − Σ)e, e = Σ̂ − Σ
sup
n e∈E i=1
e∈E

∞

,

n
n
o
o
n
n
1X
1X
2
2
≤
sup E (fi⊗2 − Σ)e1 , e1
sup E fi⊗2 − Σ, e1 ⊗ e2
n i=1 e1 ∈E
n i=1 e1 ,e2 ∈E
n



1X
sup
E fi⊗4 − Σ⊗2 (e1 ⊗ e2 ), e1 ⊗ e2 = E fi⊗4 − Σ⊗2
n i=1 e1 ,e2 ∈E

∞

.

2

∞

= kfi kL2 ≤ c2 = U , then σ 2 ≤ U 2 .

The weak variance for Gaussian data


Similarly to the bounded case, we estimate E f ⊗4 − Σ⊗2 p for Gaussian data. Consider f from a Gaussian process
with mean zero and covariance Σ and define fs = f (s). Strictly speaking these variables are not norm bounded, but similar
to concentration results for Gaussian random variables in Rd (Giné and Nickl, 2016, Theorem 3.1.9), we found below that
our methods still work well. The integral kernel can be written as (Isserlis, 1918)
E (fs ft fu fv ) = E (fs ft ) E (fu fv ) + E (fs fu ) E (ft fv ) + E (fs fv ) E (ft fu )
= cf (s, t)cf (u, v) + cf (s, u)cf (t, v) + cf (s, v)cf (t, u).

4


Hence, we have that E (fs ft fu fv ) − Σs,t Σu,v = Σs,u Σt,v + Σs,v Σt,u and that the operator E f ⊗4 − Σ⊗2 , which can be
thought of as an Hilbert-Schmidt operator on the space Op(L2 ), can be represented by the integral kernel cf (s, u)cf (t, v) +
cf (s, v)cf (t, u). These two terms are merely relabeled versions of Σ⊗2 . Consequently, using the subadditivity of the norm,
E f ⊗4 − Σ⊗2 p ≤ Σ⊗2 p + Σ⊗2 p = 2 Σ⊗2 p . For example, for the Hilbert-Schmidt norm,

E f ⊗4 − Σ⊗2

2
HS

ZZZZ

2

{cf (s, u)cf (t, v) + cf (s, v)cf (t, u)} dsdtdudv
ZZZZ
4
4
= 2 kΣkHS + 2
cf (s, u)cf (s, v)cf (t, v)cf (t, u)dsdtdudv ≤ 4 kΣkHS .

=

Lemma 5.1 of Horváth and Kokoszka (2012) gives an explicit form of a covariance operator of Σ in terms of the
eigenfunctions of Σ for Gaussian data in the Hilbert-Schmidt setting.
Given λi , the eigenvalues of Σ, the spectrum of Σ⊗2 is {λi λj }∞
i,j=1 . Hence, for any of the p-Schatten norms,
2
kΣ ⊗ Σkp = kΣkp . Note that in the above calculations, the weak variance depends on the unknown Σ. In practice, this can
be replaced by the empirical estimate Σ̂.

4

Applications

4.1 k Sample Comparison
Testing for the equality of means among multiple sets of data is a common task in data analysis. In the functional setting,
there has been recent work on performing such a test on covariance operators in order to test whether or not k sets of curves
have similar variation. Panaretos et al. (2010) propose such a method for a two sample test on covariance operators given
data from Gaussian processes. Similarly, Fremdt et al. (2013) propose a non-parametric two sample test on covariance
operators. Both of these approaches make use of the Karhunen-Loève expansion and, hence, the underlying Hilbert space
geometry. Pigoli et al. (2014) take a comparative look at a variety of metrics to rank their statistical power when used in a
two sample permutation test.
Following from the results of Pigoli et al. (2014), our method uses the p-Schatten norms with the concentration
inequality based confidence sets of the previous section to compare covariance operators. In the two sample setting, we
are able to achieve similar statistical power to that of the permutation test after proper tuning of the coefficients in the
inequalities. Furthermore, the analytic nature of the concentration approach leads to a significant reduction in computing
time, which offers an even more significant savings for larger values of k.
From the confidence set constructed in the previous section, we can devise a test for comparing the empirical covariance
(1)
(1)
(k)
(k)
operators generated from k samples of functional data. Let the k samples be f1 , . . . , fn1 , . . . , f1 , . . . , fnk where for
(i)
each sample i and all elements j = 1, . . . , ni , fj has covariance Σ(i) . Our goal is to design a test for the following two
hypotheses:
H0 : Σ(1) = . . . = Σ(k)

H1 : ∃ i, j s.t. Σ(i) 6= Σ(j) .

To achieve this, a pooled estimate of the weak variance is computed as a weighted average of each sample’s individual
weak variance in similar style to that of a standard t-test. Let the total data size be N = n1 + . . . + nk and σi2 be the
Pk
2
weak variance for sample i, then the pooled variance is defined as σpool
= N −1 i=1 ni σi2 . Given Gaussian data and the
Pk
2
2
p-Schatten norm, for example, this reduces to σpool
= 2N −1 i=1 ni kΣ(i) k2p . In practice, σpool
is estimated from the data
for the following confidence regions in order to have those regions only depend on the data.
Taking inspiration from the standard analysis of variance (Casella and Berger, 2002, Chapter 11), let Σ̂(i) be the
empirical estimate of the covariance operator for the ith sample, and let Σ̂ be the estimate of the covariance operator for the
total data set. Making use of the confidence sets for covariance operators from Section 3.2 gives the rejection region

C=





f:

k
X
i=1

Σ̂(i) − Σ̂

>
p

ni
k
X
X
i=1



(i) ⊗2
εi,j fj
− Σ̂

j=1

+
p

+

5

k
X

2
σpool

i=1

ni

!1/2
1/2

(−2 log 2α)

+

k
X
σpool
i=1

ni

!


log 2α 
3



,

which under the null hypothesis will have size no greater than the desired α.
The size of the test induced by this rejection region is significantly less than the target size α due to the use of multiple
concentration inequalities. Hence, tuning the inequalities is required to yield a useful test. Many experiments were run on
simulated data sets generated as samples from a Gaussian process with randomly generated covariance operators whose
eigenvalues were chosen to decay at a variety of rates. In this setting, the coefficients of 1 − k −1/2 for the Rademacher term
and (k + 2)/(k + 3) for the deviation term were determined experimentally to improve the size of the confidence region:

ni
k
k

X


X
X
(i) ⊗2
C = f :
Σ̂(i) − Σ̂ > 1 − k −1/2
εi,j fj
− Σ̂
+
p

i=1

i=1

+

4.2

j=1

p


!1/2

 X
k
2
σpool
k+2
k+3



i=1

ni

(−2 log 2α)1/2 +

k
X
σpool
i=1

!

ni


log 2α 
. (4.1)
3 

Classification of Operators

Classification of functional data has been an area of heavy research over the last two decades. James and Hastie (2001)
extend linear discriminant analysis to functional data. Hall et al. (2001) and Glendinning and Herbert (2003) classify with
principle components. Ferraty and Vieu (2003) implement kernel estimators. General linear models for functional data are
discussed by Müller and Stadtmüller (2005). Delaigle and Hall (2012) analyze the asymptotic properties of the centroid
based classifier. Wavelet based classification is detailed by Berlinet et al. (2008) and Chang et al. (2014).
One application of our method beyond classification of functional data is the classification of covariance operators. In
the setting of speech analysis, consider multiple speakers and multiple samples of speech from each speaker. The speech
samples can be combined into a single sample covariance operator for each speaker. Then, our method can be employed,
for example, to classify the covariance operators by speaker gender or speaker language. Evidence that this is a fruitful
approach can be found in the analysis of Pigoli et al. (2014) and Pigoli et al. (2015) where a variety of metrics are compared
for their efficacy when performing inference on covariance operators. These articles detail the discrepancy between sample
covariance operators produced by speakers of different romance languages.
Given k possible labels and n samples of labeled data (Yi , fi ) with label Yi ∈ {1, . . . , k} and observation fi ∈ L2 (I),
our goal is to determine the probability that a newly observed g ∈ L2 (I) belongs to label Y = j. Given such a g, the Bayes
classifier chooses the label y = arg maxj pr(Y = j | g) where pr(Y = j | g) = pr(g | Y = j)pr (Y = j)/pr (g).
Beginning
mean of each category iscomputed:
P with a training set of n samples with nj samples of label j, the sample
¯
¯
¯
} + r]
f¯j = n−1
j
i:Yi =j fi . The probability pr(g | Y = j) above is replaced with pr[ fj − g L2 > E{ fj − E fj
L2

¯
¯
with the goal of making a decision based on how much more fj differs from g than fj differs from its expectation E f¯j .
Similar techniques to those is Section 3 as used. Define the Rademacher sum, Rj , and the empirical weak variance, σ̂j2 , for
label j to be, respectively,
Rj =

1 X
εi (fi − f¯j ),
nj

σ̂j2 =

i:Yi =j

1 X ⊗2 ¯⊗2
fi − fj
nj
i:Yi =j

p

where εi are independent and identically distributed Rademacher random variables. The tail bound for the above probability
is then
!

−nj r2
¯
pr fj − g L2 − kRj kL2 > r < exp
,
(4.2)
4 kRj kL2 U + 2σ̂j2 + 2rU/3

where U is an upper bound on kfi kL2 . However, this can be approximated by the Gaussian tail exp −nj r2 /2σj2 . In the
simulations of Section 5.3, this approximation actually achieves a better correct classification rate on both Gaussian and
t-distributed data. This specifically works on t-distributed data as the estimate in Equation 4.3 below is merely concerned
with comparing the tail bounds rather than their specific values. Consequently, the tail for every category is underestimated
in the t case, but the ratio remains valid for comparison purposes.
Assuming uniform priors on the labels, the estimate for the probability expression in the Bayes classifier is

!2 


¯
−
kR
k
f
−
g
φj (g)
nj
j
j L2
L2
pr(Y = j | g) ≈ Pk
,
φj (g) = exp −
.
(4.3)
 2

σ̂j
φl (g)
l=1

6

This can be extended to the case where an unlabeled observation is a collection of curves g1 , . . . , gm by replacing
kf¯j − gkL2 in the above expression with kΣ̄j − Σ̂g kp where Σ̄j is the sample covariance of the fi with label j and Σ̂g is
the sample covariance of the gi . The Rademacher and weak variance terms would also be updated accordingly.

4.3

Clustering of Operator Mixtures

Closely related to the problem of classification is the problem of clustering. Given a sample of functional data, we want to
assign one of a finite collection of labels to each curve. For example, in speech processing, one may want to cluster sound
clips based the language of the speaker, or, to be discussed in Section 5.4, one may want to separate unlabeled phoneme
curves into clusters.
There have been many recently proposed methods for clustering functional data. Many approaches begin by constructing
a low dimensional representation of the data in some basis such as modelling the data with a B-spline basis followed
by clustering the spline representations with k-means (Abraham et al., 2003). A similar approach makes use of the
eigenfunctions of the covariance operator instead of B-splines (Peng and Müller, 2008). In contrast, we will attempt
to cluster functions or operators directly via a concentration of measure approach similar to the previously described
classification procedure.
Consider the same setting to the previous section of multiple observations from multiple categories. However, now
the category labels are missing. This is a functional mixture model where each observed functional datum is a stochastic
process with one of k possible covariance operators. In the below experiments, the data will be simulated from a Gaussian
process. The goal is to correctly separate the data into k sets. To achieve this, an expectation-maximization style algorithm
is implemented.
(i)
(i)
Let the observed operator data be S1 , . . . , Sn ∈ Op(L2 ) where each Si = cov(f1 , . . . , fmi ) is a rank mi operator
produced from mi functional observations. Let the latent label variables be Y1 , . . . , Yn ∈ {1, . . . , k}. Assuming no prior
knowledge on the proportions of data in each category, the algorithm is initialized with the Jeffreys prior for the Dirichlet
(0)
distribution by randomly generating ρi,· ∼ Dirichlet (1/2, . . . , 1/2) , the initial probability vector that pr(Yi = ∗ | fi ).
(t)

Assuming t iterations of the algorithm have completed, we have a label probability vector ρi,· for each of the n
(t+1)

observations. Given this collection of vectors, the expected proportions of each category can be estimated as τj
=
Pn
Pn
(t)
(t+1)
−1
−1
n
, and a weighted Rademacher sum,
i=1 E (1Yi =j ) = n
i=1 ρi,j . Similarly, a weighted sum of the data, Σ̂j
(t+1)

Rj

, can be used to update the estimated covariance operators for each label j:
(t+1)

Σ̂j

Pn

Pn

(t)

i=1 ρi,j Si
,
= Pn
(t)
i=1 ρi,j

(t+1)

Rj

=

i=1



(t)
(t+1)
ρi,j εi Si − Σ̂j
.
Pn
(t)
i=1 ρi,j

Lastly, a pooled weak variance is required, which is used in place of each individual category weak variance. Otherwise,
(t+1)
in practice, one single category captures all of the data points. By defining the pooled covariance operator as Σ̂pool =
Pk
(t+1) (t+1)
(t+1)
Σ̂j
, then the pooled weak variance in the Gaussian case, for example, is estimated by 2kΣ̂pool kp .
j=1 τj
(t)

As a result, the label probability vectors ρi,· can be updated given the t + 1st collection of estimated covariance
operators, Rademacher sums, and the pooled covariance operator. From the previous section, Equation 4.3 can be used
(t+1)
(t+1)
(t+1)
to determine ρi,j = pr(Yi = j | Si , Σ̂1
, . . . , Σ̂k
), the probability that observation i belongs to the jth category.
This process can be iterated until a local optimum is reached.

5
5.1

Numerical Experiments
Simulated and phoneme data

To test each of the above three applications, experiments were first run on simulated data. These data sets were generated as
observations from Gaussian or t-distributed processes. The covaraiance operators were randomly chosen given a specific
decay rate for the eigenvalues.
Secondly, the phoneme data to be tested (Ferraty and Vieu, 2003; Hastie et al., 1995) is a collection of 400 logperiodograms for each of five different phonemes: /A/ as in the vowel of “dark”; /O/ as in the first vowel of “water”; /d/
as in the plosive of “dark”; /i/ as in the vowel of “she”; /S/ as in the fricative of “she”. Each curve contains the first 150
frequencies from a 32 ms sound clip sampled at a rate of 16-kHz.
7

5.2

k Sample Comparison

The above confidence set in Equation 4.1 comparing k samples can be used to refute the null hypothesis that all covariance
operators are equal. A two sample permutation test was performed in Pigoli et al. (2014). Given two samples of
(1)
(1)
(2)
(2)
functional data, f1 , . . . , fn and f1 , . . . , fm with associated covariance operators Σ(1) and Σ(2) , respectively, the
desired hypotheses to test are
H0 : Σ(1) = Σ(2)

H1 : Σ(1) 6= Σ(2) .

When using a permutation test, the labels are randomly reassigned M times, and each time, the distance between the two
new covariance operators is computed. For sufficiently large M , this procedure will return the exact significance level of
the observations with respect to the data set.
A power analysis was performed between the permutation method and our proposed concentration approach using Equation 4.1. Given two different operators Σ(1) and Σ(2) and γ ≥ 0, an interpolation between the two operators is constructed
as Σ(γ) = [(Σ(1) )1/2 + γ{S(Σ(2) )1/2 − (Σ(1) )1/2 }][(Σ(1) )1/2 + γ{S(Σ(2) )1/2 − (Σ(1) )1/2 }]∗ , where S is an operator
minimizing the Procrustes distance, between Σ(1) and Σ(2) , which is dProc (Σ(1) , Σ(2) )2 = inf S∈U {L2 (I)} kR(1) − R(2) Sk22

where Σ(i) = (R(i) )(R(i) )∗ and U L2 (I) is the space of unitary operators on L2 (I) (Pigoli et al., 2014).
Monte Carlo simulations were run in order to estimate the power of each test. Two operators Σ(1) and Σ(2) with
similar eigenvalue decay were compared. For sample size n = 50, γ ∈ {0, .1, .2, .3, .4, .5}, and for sample size n = 500,
γ ∈ {0, .05, .1, .15, .2, .25}. For each γ, 5000 samples of size n were generated for Σ(1) and Σ(γ) . Equation 4.1 and the
permutation method (Pigoli et al., 2014) were both implemented to estimate the empirical power.
Figures 1 and 2 display the results for operators whose eigenvalues decay at a quartic and quadratic rate, respectively.
The solid lines indicate the power of the permutation test, and the circle lines indicate the power of our concentration
approach. The colors red, green, and blue correspond to the three norms trace, Hilbert-Schmidt, and operator, respectively.
In most cases, the concentration approach is able to achieve the same power to reject the null as does the permutation
test. The notable exception is for the trace norm when the eigenvalues decay slowly. The added benefit to the concentration
approach is the speed with which it executes. Across all of the Monte Carlo simulations, our concentration approach ran
on average 28.14 times faster than the permutation method. This was computed by tracking the amount of computation
time each method spent while producing the plots in Figures 1 and 2, which corresponds to 5 values of γ, 2 values of α, 3
different norms, 2 sample sizes, and 5000 replications each resulting in 300,000 function calls for both the permutation and
concentration methods. Unlike the other norms, the Hilbert-Schmidt norm can be calculated without explicit computation
of the eigenvalues. For each evaluation of the permutation test, 100 permutations of the data were generated, which
corresponds to 100 random draws and 100 eigenvalue computations. More accuracy would require even more permutations.
In comparison, our concentration approach requires only 3k eigenvalue computations and no random draws and hence is
only dependent on the number of samples regardless of data size or α.
The proposed k-sample test was also used to compare samples of log-periodogram curves from the spoken phonemes
/A/ and /O/ . As one can imagine, these vowels can be hard to distinguish; see Section 5.4 for further evidence of this.
For k ∈ {2, 3, 4, 5, 6}, k − 1 disjoint sets of 40 /A/ curves and one set of 40 /O/ curves were randomly sampled from the
data set. This was replicated 500 times, and each time Equation 4.1 was used to decide whether or not the k covariance
operators were equivalent at the α = 0.05 level. The resulting estimated statistical power for each k is
k
Power

2
0.00

3
0.018

4
0.228

5
0.656

6
0.936

In the null setting, the above experiement was rerun except that every disjoint set of curves came from the /A/ set. The
resulting experimentally computed test sizes are
k
Size

5.3

2
0.00

3
0.00

4
0.00

5
0.004

6
0.072

Binary and trinary classification

Our concentration of measure (CoM) method is implemented on covariance
Pn operators making use of the trace norm k·ktr
where for a covariance operator Σ with eigenvalues {λi }∞
i=1 , kΣktr =
i=1 |λi |. The trace norm was chosen based on
the analysis of the preceding section as well as that of Pigoli et al. (2014) where it achieved the best performance when
compared with the other p-Schatten norms. The CoM approach to classification of operators is tested in a variety of
simulations against other standard approaches to functional classification. The methods used for comparison are k-nearest
8

Sample Size: 500
●
●

1.0

1.0

Sample Size: 50
●
●

●
●

●

●

●

●

●

0.15

0.20

0.25

●
●

●

●

●

●

0.15

0.20

0.25

●

●

●

0.8

0.8

●

●

0.6

●
●

0.4

Empirical Power

0.6

●

0.4

Empirical Power

●

0.2

0.2

●
●
●
●

●
●
●
●

●

0.0

0.1

0.2

●

●

●

0.3

0.4

0.5

0.0

●

●

●

●

0.00

0.05

0.10

Gamma

Sample Size: 50

Sample Size: 500
●

1.0

Gamma

1.0

0.0

●
●

●

●

●

●

0.8

0.8

●

●

●

0.6

●

0.4

0.6

Empirical Power

●

0.4

Empirical Power

●

●
●

0.2

0.2

●

●
●
●

●
●
●

●

●

0.0

0.1

0.2

●

●

●

0.3

0.4

0.5

0.0

0.0

●
●

Gamma

●
●

●

●

0.00

0.05

0.10

Gamma

Figure 1: A power analysis for testing whether or not operator Σ(1) = Σ(γ) comparing the permutation method (solid
lines) with the concentration approach (circle lines). The size α = 0.05 in the top row, and α = 0.01 in the bottom. The
eigenvalues of the operators decay at a rate O(k −4 ). The red, green, and blue lines respectively correspond to the trace
class, Hilbert-Schmidt, and operator norms.

9

Sample Size: 500
●

1.0

1.0

Sample Size: 50
●
●

●
●
●

●

●

●

●

0.15

0.20

0.25

●

●
●

●

●

●
●

0.8

0.8

●

●

●

0.6

●

●

0.4

Empirical Power

0.6

●

0.4

Empirical Power

●

●

0.2

0.2

●
●

●
●
●

●

●

0.0

0.1

0.2

●

●

●

0.3

0.4

0.5

0.0

●
●
●
●

●
●
●
●

●

●

0.00

0.05

0.10

Gamma

Sample Size: 50

Sample Size: 500
1.0

Gamma

1.0

0.0

●
●
●

●

●

●

●

●
●

0.8

0.8

●

0.6

●

0.4

Empirical Power

0.6

●

0.4

Empirical Power

●
●

●

●
●

0.2

0.2

●

●
●
●
●

●

0.0

0.1

0.2

●
●
●

●

●

0.3

0.4

0.5

0.0

0.0

●
●
●
●
●

Gamma

●
●

●
●

●

0.00

0.05

0.10

●

●

●

0.15

0.20

0.25

Gamma

Figure 2: A power analysis for testing whether or not operator Σ(1) = Σ(γ) comparing the permutation method (solid
lines) with the concentration approach (circle lines). The size α = 0.05 in the top row, and α = 0.01 in the bottom. The
eigenvalues of the operators decay at a rate O(k −2 ). The red, green, and blue lines respectively correspond to the trace
class, Hilbert-Schmidt, and operator norms.

10

Table 1: A comparison of the performances of five classification methods. The first entry for each method corresponds to
classifying the covariance operators as functions. The prime entry corresponds to classifying curves with a majority vote.
The estimated percent of correct classification is listed in the table with the sample standard deviation in brackets. The
top block comes from Gaussian process data, and the bottom comes from t-process data with 4 degrees of freedom. The
highest percentage of each column is marked in bold.
k
CoM
KNN
KNN0
Kernel
Kernel0
GLM
GLM0
Tree
Tree0

1
0.62 (0.05)
0.52 (0.04)
·
0.65 (0.04)
·
0.51 (0.04)
·
0.57 (0.04)
·

2
0.62 (0.05)
0.44 (0.04)
0.47 (0.03)
0.64 (0.05)
0.70 (0.04)
0.62 (0.04)
0.50 (0.04)
0.54 (0.04)
0.55 (0.04)

4
0.76 (0.08)
0.57 (0.04)
0.59 (0.04)
0.75 (0.03)
0.82 (0.02)
0.74 (0.04)
0.50 (0.03)
0.59 (0.03)
0.59 (0.04)

8
0.87 (0.06)
0.76 (0.04)
0.74 (0.03)
0.87 (0.03)
0.92 (0.02)
0.86 (0.02)
0.50 (0.04)
0.66 (0.04)
0.60 (0.05)

16
0.96 (0.03)
0.91 (0.02)
0.89 (0.02)
0.96 (0.02)
0.99 (0.01)
0.94 (0.02)
0.50 (0.04)
0.75 (0.03)
0.60 (0.09)

CoM
KNN
KNN0
Kernel
Kernel0
GLM
GLM0
Tree
Tree0

0.59 (0.05)
0.42 (0.04)
·
0.65 (0.04)
·
0.50 (0.03)
·
0.54 (0.04)
·

0.62 (0.05)
0.45 (0.04)
0.45 (0.04)
0.64 (0.05)
0.68 (0.04)
0.62 (0.03)
0.50 (0.03)
0.54 (0.04)
0.54 (0.04)

0.75 (0.07)
0.58 (0.04)
0.58 (0.04)
0.75 (0.04)
0.80 (0.03)
0.74 (0.04)
0.50 (0.03)
0.59 (0.03)
0.57 (0.04)

0.86 (0.06)
0.76 (0.03)
0.72 (0.03)
0.87 (0.02)
0.92 (0.02)
0.86 (0.03)
0.51 (0.04)
0.67 (0.03)
0.59 (0.05)

0.95 (0.04)
0.92 (0.02)
0.89 (0.03)
0.96 (0.02)
0.99 (0.01)
0.94 (0.02)
0.50 (0.03)
0.75 (0.03)
0.57 (0.06)

CoM, concentration of measure; KNN, k-nearest-neighbours; GLM, generalized linear model; Tree, regression tree.
neighbours (Ferraty and Vieu, 2006), classification using kernel estimators (Ferraty and Vieu, 2003), general linear model
(Müller and Stadtmüller, 2005), and regression trees.
The first simulation asks each method to classify observed mean zero Gaussian process data or mean zero t-process
data with 4 degrees of freedom. The two covariance operators in question, Σ1 and Σ2 , are the sample covariances of the
male and of the females of the Berkeley growth curve data (Ramsay and Silverman, 2005). In particular, n collections of k
curves were generated from each of Σ1 and Σ2 as a training set, and m collections of k curves were generated as a test
set. The CoM method was trained on the set of n sample covariances and used to classify each of the m test covariances.
The remaining classification methods were trained and tested in two separate ways: By treating each sample covariance
as a function and classifying as usual, and by training on all n × k observations and testing each of the m collections by
classifying each constituent curve individually and taking a majority vote with ties settled by a uniform random draw.
For group sizes k = 1, 2, 4, 8, 16, 100 simulations were run with n = 100 sets of k training curves. To compare the
accuracy of each approach m = 100 sets of k testing curves were generated for each operator. The accuracy of each method
is tabulated in Table 1.
The concentration method performed well against the alternatives. Its performance was on par with the kernel method
applied to each covariance operator as a function. Our method was only consistently outperformed by the kernel method
implementing the majority vote approach. However, the two operators in question have very similar weak variances. The
next simulation demonstrates how the concentration method adapts naturally when the variances of each label significantly
differ.
Continuing from the previous simulation, a third operator is constructed from Σ1 and Σ2 by averaging these two and
then scaling up the non-principle eigenvalues by a factor of 5. This, in some sense, creates a third operator between the first
two, but with higher variance. The simulation is carried out precisely as before, but incorporating all three operators. In this
setting, our concentration approach demonstrates the best performance. The results are listed in Table 2.
These five methods tested on simulated data were also tested against phoneme data. Across 50 iterations, each set
of 400 curves was partitioned at random into an 100 curve training set and a 300 curve testing set. The five classifiers
were trained and run on each of the 300 × 5 curves individually. For our concentration of measure approach, the rank
one operator associated to each individual curve was compared with the covariance operator formed from the 100 × 5

11

Table 2: A comparison of the performances of five classification methods as in Table 1, but with three potential classes from
which to choose. The estimated percent of correct classification is listed in the table with the sample standard deviation in
brackets. The top block comes from Gaussian process data, and the bottom comes from t-process data with 4 degrees of
freedom. The highest percentage of each column is marked in bold.
k
CoM
KNN
KNN0
Kernel
Kernel0
GLM
GLM0
Tree
Tree0

1
0.51 (0.04)
0.50 (0.03)
·
0.54 (0.03)
·
0.36 (0.04)
·
0.44 (0.03)
·

2
0.55 (0.04)
0.52 (0.03)
0.55 (0.03)
0.52 (0.03)
0.58 (0.03)
0.41 (0.04)
0.35 (0.04)
0.44 (0.03)
0.46 (0.03)

4
0.75 (0.05)
0.61 (0.03)
0.68 (0.03)
0.64 (0.03)
0.69 (0.02)
0.49 (0.04)
0.36 (0.04)
0.45 (0.03)
0.51 (0.04)

8
0.89 (0.05)
0.75 (0.03)
0.80 (0.02)
0.77 (0.03)
0.81 (0.03)
0.57 (0.03)
0.36 (0.05)
0.50 (0.03)
0.51 (0.07)

16
0.97 (0.03)
0.90 (0.02)
0.90 (0.02)
0.92 (0.02)
0.92 (0.02)
0.65 (0.03)
0.35 (0.05)
0.55 (0.03)
0.47 (0.07)

CoM
KNN
KNN0
Kernel
Kernel0
GLM
GLM0
Tree
Tree0

0.46 (0.05)
0.46 (0.03)
·
0.50 (0.03)
·
0.35 (0.03)
·
0.42 (0.03)
·

0.50 (0.06)
0.49 (0.03)
0.50 (0.03)
0.48 (0.03)
0.53 (0.03)
0.41 (0.04)
0.35 (0.03)
0.44 (0.03)
0.43 (0.03)

0.63 (0.05)
0.57 (0.03)
0.64 (0.03)
0.57 (0.03)
0.66 (0.03)
0.46 (0.04)
0.36 (0.04)
0.45 (0.03)
0.47 (0.04)

0.75 (0.08)
0.67 (0.03)
0.77 (0.02)
0.68 (0.03)
0.77 (0.02)
0.53 (0.03)
0.36 (0.04)
0.46 (0.03)
0.48 (0.06)

0.85 (0.07)
0.78 (0.02)
0.87 (0.02)
0.80 (0.02)
0.85 (0.02)
0.58 (0.04)
0.36 (0.06)
0.49 (0.03)
0.46 (0.08)

CoM, concentration of measure; KNN, k-nearest-neighbours; GLM, generalized linear model; Tree, regression tree.
Table 3: Percentage of correct classification of the five phonemes against the five methods. The highest percentage of each
column is marked in bold.

CoM
KNN
Kernel
GLM
Tree

/A/
0.769
0.724
0.720
0.790
0.708

/O/
0.768
0.791
0.805
0.723
0.694

/d/
0.966
0.985
0.984
0.982
0.956

/i/
0.985
0.974
0.972
0.959
0.878

/S/
0.994
1.000
0.999
0.992
0.926

training curves. The results are detailed in Table 3. Our concentration of measure approach only uniformly outperforms the
regression tree classifier, but has comparable performance to the other three methods, and none of the competing methods
uniformly outperforms ours.

5.4

The expectation-maximization algorithm in practice

The experiments described and depicted below make use of the trace norm only. It was determined through experimentation
that the expectation-maximization algorithm we propose in Section 4.3 does not perform well under the topology of either
the Hilbert-Schmidt or operator norms as they give more emphasis to the principle eigenvalue at the expense of the others.
The usual behavior under these norms is for all estimates to converge to the average of all of the data points. This is in
contrast to the better performance of the algorithm making use of the trace norm, which is somewhat more uniform in its
treatment of the eigenstructure.
As a first test case, this algorithm was run given three target covariance operators, which were constructed by taking
three randomly generated orthonormal bases Ui and a diagonal operator D of eigenvalues decaying at a rate λk = O(k −4 )
and multiplying Σi = Ui DUi T . Let the three target covariance operators be denoted as Σa , Σb , and Σc . For each of these
operators, 1000 rank four data points were generated from a zero mean Gaussian process. From the data, the algorithm
(t)
(t)
(t)
initializes three estimates Σ̂1 , Σ̂2 , and Σ̂3 , which attempt to locate the three target operators as the method iterates.

12

Table 4: Clustering of simulated operators.

Label a
Label b
Label c

Rank 4 Operators
Cluster 1 Cluster 2 Cluster 3
0
0
1000
0
1000
0
1000
0
0

Rank 1 Operators
Cluster 1 Cluster 2 Cluster 3
169
557
274
184
549
267
483
230
287

Table 5: Clustering 2000 phoneme curves into 5 clusters.

/A/
/O/
/d/
/i/
/S/

Cluster 1
3
0
0
212
4

Cluster 2
0
0
0
0
393

Cluster 3
397
397
0
0
0

Cluster 4
0
1
378
1
0

Cluster 5
0
2
22
187
3

After 15 iterations, the original 3000 data points were perfectly separated into three groups. To make the problem harder, a
second test case was run identical to the first except that the observed operators are all of rank one. The resulting clusters
from both tests are in Table 4. The inaccuracy in the rank one setting is equivalent to the poor performance of classification
of rank one operators detailed in Tables 1 and 2.
For the phoneme data all 400 sample curves from each of the five phoneme sets were clustered individually as curves.
The algorithm was run for 20 iterations and told to partition the data into five clusters. The results are in Table 5. Cluster 3
captured almost all of the vowels /A/ and /O/ , which, recalling their definition in Section 5.1, are quite similar in sound.
Clusters 2 and 4 contain the majority of /S/ and /d/ curves, respectively. Lastly, Clusters 1 and 5 split the set of /i/
curves roughly in half. Rerunning the experiment with only four clustered nicely separated the data as displayed in Table 6.
Cluster 3 again contains almost all of the /A/ and /O/ . The remaining clusters mostly partition the /d/ , /i/ , and
/S/ with high accuracy. Hence, barring the ability to distinguish the very similar /A/ and /O/ curves, the proposed
expectation-maximization algorithm is an effective method for the unsupervised clustering of phonemes.

References
Christophe Abraham, Pierre-André Cornillon, ERIC Matzner-Løber, and Nicolas Molinari. Unsupervised curve clustering
using b-splines. Scandinavian journal of statistics, 30(3):581–595, 2003.
Sylvain Arlot, Gilles Blanchard, and Etienne Roquain. Some nonasymptotic results on resampling in high dimension, i:
confidence regions. The Annals of Statistics, 38(1):51–82, 2010.
Peter L Bartlett and Shahar Mendelson. Rademacher and gaussian complexities: Risk bounds and structural results. The
Journal of Machine Learning Research, 3:463–482, 2003.
Peter L Bartlett, Stéphane Boucheron, and Gábor Lugosi. Model selection and error estimation. Machine Learning, 48
(1-3):85–113, 2002.
Alain Berlinet, Gérard Biau, and Laurent Rouviere. Functional supervised classification with wavelets. In Annales de
l’ISUP, volume 52, 2008.
Table 6: Clustering 2000 phoneme curves into four clusters.

/A/
/O/
/d/
/i/
/S/

Cluster 1
0
1
384
1
0

Cluster 2
0
0
1
5
397

13

Cluster 3
400
398
0
0
0

Cluster 4
0
1
15
394
3

Denis Bosq. Linear processes in function spaces: theory and applications, volume 149. Springer Science & Business
Media, 2012.
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities: A nonasymptotic theory of independence. Oxford University Press, 2013.
George Casella and Roger L Berger. Statistical inference, volume 2. Duxbury Pacific Grove, CA, 2002.
Chung Chang, Yakuan Chen, and R Todd Ogden. Functional data classification: a wavelet approach. Computational
Statistics, 29(6):1497–1513, 2014.
Aurore Delaigle and Peter Hall. Achieving near perfect classification for functional data. Journal of the Royal Statistical
Society: Series B (Statistical Methodology), 74(2):267–286, 2012.
Zhou Fan. Confidence regions for infinite-dimensional statistical parameters. Part III essay in Mathematics, University of
Cambridge, 2011. http://web.stanford.edu/˜zhoufan/PartIIIEssay.pdf.
Frédéric Ferraty and Philippe Vieu. Curves discrimination: a nonparametric functional approach. Computational Statistics
& Data Analysis, 44(1):161–173, 2003.
Frédéric Ferraty and Philippe Vieu. Nonparametric functional data analysis: theory and practice. Springer Science &
Business Media, 2006.
Stefan Fremdt, Josef G Steinebach, Lajos Horváth, and Piotr Kokoszka. Testing the equality of covariance operators in
functional samples. Scandinavian Journal of Statistics, 40(1):138–152, 2013.
Evarist Giné and Richard Nickl. Confidence bands in density estimation. The Annals of Statistics, 38(2):1122–1170, 2010.
Evarist Giné and Richard Nickl. Mathematical Foundations of Infinite-Dimensional Statistical Models. Cambridge
University Press, 2016.
Richard H Glendinning and RA Herbert. Shape classification using smooth principal components. Pattern recognition
letters, 24(12):2021–2030, 2003.
Peter Hall, DS Poskitt, and Brett Presnell. A functional dataanalytic approach to signal discrimination. Technometrics, 43
(1):1–9, 2001.
Trevor Hastie, Andreas Buja, and Robert Tibshirani. Penalized discriminant analysis. The Annals of Statistics, pages
73–102, 1995.
Lajos Horváth and Piotr Kokoszka. Inference for functional data with applications, volume 200. Springer Science &
Business Media, 2012.
Leon Isserlis. On a formula for the product-moment coefficient of any order of a normal frequency distribution in any
number of variables. Biometrika, pages 134–139, 1918.
Gareth M James and Trevor J Hastie. Functional linear discriminant analysis for irregularly sampled curves. Journal of the
Royal Statistical Society. Series B, Statistical Methodology, pages 533–550, 2001.
Ci-Ren Jiang, John AD Aston, and Jane-Ling Wang. A functional approach to deconvolve dynamic neuroimaging data.
arXiv preprint arXiv:1411.2051, 2014.
Gerard Kerkyacharian, Richard Nickl, and Dominique Picard. Concentration inequalities and confidence bands for needlet
density estimators on compact homogeneous manifolds. Probability Theory and Related Fields, 153(1-2):363–404,
2012.
Vladimir Koltchinskii. Rademacher penalties and structural risk minimization. Information Theory, IEEE Transactions on,
47(5):1902–1914, 2001.
Vladimir Koltchinskii. Local rademacher complexities and oracle inequalities in risk minimization. The Annals of Statistics,
34(6):2593–2656, 2006.
Michel Ledoux. The concentration of measure phenomenon, volume 89. American Mathematical Soc., 2001.
14

Karim Lounici and Richard Nickl. Global uniform risk bounds for wavelet deconvolution estimators. The Annals of
Statistics, 39(1):201–231, 2011.
Hans-Georg Müller and Ulrich Stadtmüller. Generalized functional linear models. Annals of Statistics, pages 774–805,
2005.
Victor M Panaretos, David Kraus, and John H Maddocks. Second-order comparison of gaussian random functions and the
geometry of dna minicircles. Journal of the American Statistical Association, 105(490):670–682, 2010.
Jie Peng and Hans-Georg Müller. Distance-based clustering of sparsely observed stochastic processes, with applications to
online auctions. The Annals of Applied Statistics, pages 1056–1077, 2008.
Davide Pigoli, John AD Aston, Ian L Dryden, and Piercesare Secchi. Distances and inference for covariance operators.
Biometrika, page asu008, 2014.
Davide Pigoli, Pantelis Z Hadjipantelis, John S Coleman, and John AD Aston. The analysis of acoustic phonetic data:
exploring differences in the spoken romance languages. arXiv preprint arXiv:1507.07587, 2015.
James O Ramsay and Bernard W Silverman. Functional data analysis. New York: Springer, 2005.
Walter Rudin. Real and complex analysis. Tata McGraw-Hill Education, 1987.
Michel Talagrand. New concentration inequalities in product spaces. Inventiones mathematicae, 126(3):505–563, 1996.

15

