On Perfect Clustering for Gaussian Processes
Authors: Juan Cuesta-Albertos
Subhajit Dutta
Venue: Accepted by TMLR ()

Abstract: In this paper, we propose a data based transformation for infinite-dimensional Gaussian processes and derive its limit theorem. For a clustering problem using mixture models, an appropriate modification of this transformation asymptotically leads to perfect separation of the populations under rather general conditions, except the scenario in which differences between clusters depend only on the locations; in which case our procedure is useless. Theoretical properties related to label consistency are studied for the k-means clustering algorithm when used on this transformed data. Good empirical performance of the proposed methodology is demonstrated using simulated as well as benchmark data sets, when compared with some popular parametric and nonparametric methods for such functional data.

Published in Transactions on Machine Learning Research (09/2023)
On Perfect Clustering for Gaussian Processes
Juan A. Cuesta-Albertos Departamento de Matemáticos, Estadística Computación Universidad de Cantabria, Spain
cuestajShumionn.cs
Subhajit Dutta
duttas @hitk. ac in
Deportment Mathematics ond Statistics IIT Kampur, India
Reviewed on OpenReview: https: // openreview. net/ forum? id- 4gDOV2KBuM
Abstract
In this paper, we propose a data based transformation for infinite-dimensionaal Gaussian processes and derive its limit theorem. For a clustering problem using mixture models, an appropriate modification of this transformation asymptotically leads perfect separation of the populations under rather general conditions, except the scenario in which differ- ences between clusters depend only the locations; in which case our procodure is useless. Theoretical properties ies related to label consistency are studied for the k-means clustering al- gorithm when used on this transformed data. Good empirical performance of the proposed hethodolggy is demonstrated using simulated as well as benchmark data sets, when com- pared with some popular parametric and nonparametric methods for such functional data.
1 Introduction
Suppose that we given two Gaussian distributions (say, GDs) Pi and IP2. The Hajek and Feldman property (established independently by Hajek (1958) and Feldman (1958)) states that Pi and I2 are either equivalent, else mutually singular. In other words, for every measurable set A, (A) 0 if and only if (22(A) = 0, or else there exist two disjoint measurable sets Si and S2 such that
Pi(Si) = 1, P2(()) 0 and Pi = 0, 0. = 1.
Mutual singularity is not very interesting in finite dimensions because happens only when least one of the covariance matrices is singular. However, in the functional setting, this singularity appears in non-trivial situations. To mention an example, it was shown by Rao and Varadarajan (1963) that. if the covariance operators of IP and I2, namely, y and P satisfy P = a for some a * 1, then P, and IP are mutually singular.
Tt is clear that the mutually singular case of the Hajek and Feldman property (say, HFp) looks very promising for classification as well as clustcring (in the mixture setting) of data points. Recently, some results have appeared laking advantage of this property propose perfect classifiers (see the references given below). However, the clustering problem seems to be harder and far available paper with results in this area. The main drawback of the paper by Delaigle et al (2019) is that it primarily deals with location problems (see Section M of Appendix Il for a detailed discussion). In this paper, we present a family of Iransformations on functional data which allows one to identify some mutually singular situations. The transformed data are then used to obtain perfect clustering in the mixture setting.
know, Delaigle et al (2019) is the only
To give an overview of our main contributions, let us consider Gaussian process (say, Z) defined on a bounded real interval, which without loss of generality, identify with the unit interval (0, 1]. Further,

---

Published in Transactions on Machine Learning Research (09/2023)
assume that its trajectories belong to the Hilbert space of square integrable functions H. which is defined as follows:
HI set of real functions f(t) with E [o, 1] such that J. P(t)dt < co.
The inner product in HI is (..) = S. [())())d. The keystone of this paper is Theorem 2.1 in Section 2. It states that under appropriate assumptions, if b H., then the limit of sequence of scaled Mahalanobis distances between some finitedimeensioall projections of Z and b converges in probability to a non-random limit. Scaling is done using the dimension of the projection, and this convergence holds as the dimension goes to infinity, Practical interest of this result lies in the fact that the limit depends only on the distribution of Z (say, Pz). Therefore, Thcorem 2.1 allows one to identify some cases in which GDs are mutually singular. In such scenarios, this result allows one to obtain perfect classification (also see Cuesta-Alberton and Dutta (2022)) as well as perfect clustering. We now explain this point a bit more precisely.
Consider a probability distribution I such that P = EK ThPn, where 0 < Tr < 1 with Σk= Ty = 1 and PA are GDs on H for S h S J. Additionally, assume P to be known, but the precise values of J, Th and PA for ISh S J are unknoun. Thus, Pz can be understood as a two-step procedure. In the first step, we randomly select one among the J possibilities with probabilities Th for S h S J. Now, if Pho is selected in this step, then one chooses a random function using the probability Pho, for ho E {I.....J}: We will stick to this interpretation throughout the paper.
According to this mixture model, every function Z ~ P is in fact generated from one of the PA 's Using Theorem 2.1, if we project the available functions on certain subspaces of increasing dimensions, then the proposed transformation leads to different limits (depending just on the PA which generated each function). Thus, if we have a set of observations (with at least one observation from each PA), in the limit we can identify the value of J as well as the subset of the observations produced by each P2 without any possibility of mistake. We call this fact perfect clustering, which in practice means that we can estimate J and the cluster assignment of all the observations correctly with a high probability when the dimension is large.
We believe that HFp should have attracted the attention of researchers in classification and clustering for functional data, the 'orthogonality case' apparently being more attractive because would allow one to obtain perfect classification and perfect clustering. It took 50 years before the HFp was formally used in classification. To the best of our knowledge, the first paper using HFp in classification was Baillo et al (2011), where the authors derived classification procedure using likelihood ratios. They focused on the 'oquivalence case' and hence, did not obtain perfect classification. Optimal classification of Gaussian processes (say, GPs) was analyzed in Torrecilla et al (2020) from the HFP viewpoint. Further, the optimal (Bayes') classifier of equivalent GPs was derived and a procedure to obtain asymptotically perfect classification of mutually singular GPs was described as well. The results covered both homoscodastic as well as heteroscedastic cases. Additionally, Delaigle and Hall (2012) and Delaigle and Hall (2013) investigated conditions under which a perfect classification procedure for GPs was possible and developed related classifiers. The paper by Dai et al (2017) proposed a functional classifier based on ratio of density functions, which also leads to perfect classification. These papers contain reference of the HFp. In fact, the relationship between Delaigle and Hall (2012) and the HFp was analyzed in Berrendero et al (2018), where the authors presented an expression of the optimal Bayes' rule in some classification problems for functional data. As mentioned earlier, perfect clustering has been studied only by Delaigle et al (2019).
In Rao and Varadarajan (1963) and Shepp (1966b), the authors obtained characterizations of the singularity or equivalence of Gaussian measures in functional spaces. Their results involve increasing sequences of subspaces. For equivalent GDs, the limit obtained in Rao and Varadarajan (1963) includes a term which is the exponential of expression involving the difference of the means of IP and I22. Curiously, the logarithm of this term is related with the expressions of our limits. Similarities between our proposal and those in Rao and Varadarajan (1963) and Shepp (1966b) end here because the other involved terms are different. Moreover, we handle Mahalanobis distances between data points, while these papers use Hellinger and Jeffreys functionals to measure discrepancy between distributions. As a consoquence, the characterizations they obtain are not applicable in practice to cluster data points because they depend on the underlying distribution. Moreover, it is not straight forward to compute such functionals using data points.
2

---

Published in Transactions on Machine Learning Research (09/2023)
1.1 Contributions
In this paper, we first analyze the limit of the above mentioned scaled Mahalanobis' distances by assuming the underlying parameters of the GPs to be known (see Section 2). We begin with general concentration result in Theorem 2.1. Then, we propose a transformation for clustering that asymptotically yield perfect seperation among the clusters (see Theorem 2.5). Further, this transformation can be used to find the unknown number of clusters (see Proposition 2.3). In Section 3, we estimate the covariance operator of the mixture distribution from data and state related asymptotic results for the proposed transformation. In Theorem 3.2, we prove uniform (on the sample points) consistency of the empirical version for the transformation associated with GP clustering. It is surprising that our GP clustering method fails to discriminate 'location only' scenarios, but yields perfect clustering for 'differences in scales' (sce Remark 2.2.2). We have also compared our work both theoretically (see Section M of Appendix II) as well as numerically (see Sections 4, 5 and Section N in Appendix II) with the existing literature on perfect clustering for functional data. All proofs are deferred to Appendix Some additional material is presented in Appendix I1, which includes a possible extension to non- non-Gaussian distributions (sce Scction K), discussion on a clustering procodure in the 'location only' case (see Section L) and theoretical comparisons of our results with those obtained in the paper by Delaigle et al (2019) (see Section M).
In summary, the main contributions of our paper are as follows:
1. Given a sample from a mixture GD, we develop a method which allows us to identify the number of members in the mixture as well as the cluster membership of the functions.
2. On the theoretical side, if know the probability IP (but neither nor the PAS) and we have a finite sample taken from P, then the procedure allows to identify J as well as the functions produced by each PA for 1 < h <
3. From the practical point of view, we can estimate the covariance operator of the full mixture if the sample size increases. Under appropriate conditions, our procedure allows to identify J as well as the functions produced by each Pr, with a high probability for s S J.
4. The proposed procedure only works in those cases in which the distributions in the sample have different covariance operators and we have enough difference in the values between L and L for 1 < h k < J (see Theorem 2.2 and Remark 2.5.2).
5. Our simulations as well as analysis of real data sets show promising behaviour of the procodure when compared with several existing methods.
1.2 Notation
We will use the following notation. The distribution of the random process Z will be denoted as F'z, its mean function by n² and its covariance operator (referred to simply covariance) by Σ² We use E² to denote the covariance between Z(s) and Z(t) for 8, E [0, 1]. Given a square matrix A, trace(A) denotes its trace. The usual Euclidean norm on R is denoted by || I.. To simplify notation, we do not always explicitly state the dependence of the norm on the dimension d. Further, we will assume that all the involved random quantitics are defined on a common probability space (2, The notions convergence in probability and equality in distribution are denotod by & and ~ respectively.
2 Transformation with Known Gaussian Distributions
Let {VaJaeN be an increasing sequence of subspaces of H... Here, the dimension of Va is d. This restriction is not necessary for the development which follows as long as the dimension Va goes to infinity with increasing d, but it simplifies the notation. Given the subspace Va. let 2 and Σ represent the d-dimensional mean and the d covariance matrix of the projection of Z on Vg. If u E H. we denote us (ui...., ua)" to be its projection on Va-
3

---

Published in Transactions on Machine Learning Research (09/2023)
Fix b H. Theorem 2.1 analyses the behaviour of the limit of squared Mahalanobis norm of the d-dimensional random vector (Z ba)a for d N. For every positive definite (p.d.) d x d matrix Ad, we define the map
(u, v) A / v)u for H.
(1)
Later, we will use this function alongwith the sequence of covariance matrices (>2 1
In this section, the underlying distributions are assumed to be known. After stating Theorem 2.1 and some remarks related it, we will look into an application cluster analysis inspired from this result. We will take advantage of the fact that the limit in this theorem not random, but may depend on the underlying probability distribution Pz.
Theorem 2.1 Let {AaJueN be sequence of d d symmetric, p.d. matrices and a4 ................................ a be the eigenvalues of the matrir Sa = (Aa)-1/252 (Aa) for d E N. We define aa = (af. ......... ag) and ||a||xo- max(aa)........ a4) is the supremum norm. Let us assume that
|
lim d
(2)
Let b E Hi such that there erist constants Lu and Ls (finite, or not) with
Lp = lim D (n², b) and d
(3)
Ls = lim [lmee((a)) d d
(4)
Then, DA" P L: Lu + Ls as + oo.
Remark 2.1.1 condition Theorem 2.1 is required to ensure that no single component. is extremely influential. For instance, it may happen that we take a sequence such that af d and a o(d for every 2 2 < d. Under this condition, limiting value is not feasible in Thcorem 2.1. However, this possibility is excluded by assumption (2).
Remark 2.1.2 We allow both the constants in Theorem 2.1 to be infinite. When Ls is finite, Lemma B.1 (see Appendix T) shows that assumption (2) follows from assumption (4).
Remark 2.1.3 Let Zi and Z₂ be independent observations generated from the GDs PI and Pz. Thus, Z, - Zn is GP with mean 2 -p²² and covariance 21 + Σ²2, Consider the matrix S = (A,)-1/2(2) + (Aa) -1/2 with d E N. Define Z Z, Ze and b 0 in Theorem 2.1. Then, under the assumptions of this result, the following convergence result holds:
DA Z2) D Z2, 0) F. L Ly + Lg d oc.
Here, L, limg D (n21 (02,,,,2))) 22 and Ls is defined in (4) of Theorem 2.1.
Remark 2.1.4 In general, the fact that Va C Va+1 does not guarantee the existence of any relationship be- tween the sets faf....,af and falt1 ..... a 441 holds (for instance, when Va is generated by the first d eigenfunctions of 52 in Section 2.1).
d+1) However, in some cases {a (af..... ag] } C {a ..... ............. d++1
0+11
2.1 Application: Cluster Analysis
In this subsection, we deal with random function whose distribution a two component mixture distribution of the form: Pz xiPi + TzP2, where 0 < mı < 1 and n2 1. Here, P,, denotes the GD III with Z₁ Fn having mean function z and covariance Ea for h = 1,2.
The mean function and the covariance of the mixture satisfy 2z (t) = 1 (t) + Z2) with t E [0, 1] and
(s, क Z. (s, 2.² कi 22222 (s) و (s)] 2, (1) 2e (1) for 8, [0. (5)
4

---

Published in Transactions on Machine Learning Research (09/2023)
Given N independent and identically distributed (i..dd)) r.vls Zi...., ZN Pz: consider the following set: Cn = {j: Z, was generated from Py for S S N}
(6)
with h E {1, Clearly, the set Cn depends on the sample size N. The components of the mixture distribution Pz and the sets C. for h = 1,2 are unknown, and the problem we are addressing here is the estimation of these sets. However, we assume Pz and the sets C₁ and C2 to be known in this section to build the fundamental idea behind using the proposed transformation for clustering of GPs.
Let Va with d N denote the sequence of d-dimensional subspaces generated by the d eigenfunctions associated with the d largest eigenvalues of Σ² (recall the discussion in Remark 2.1.4). Consider sub- sample of size two from Zi.... ZN generated from P. To simplify our notation, we denote them to be Z and Z2, respectively. They are independent, with Izz, = Ph and 12z, = Ps for h, k E [1,2]. The clustering procedure that propose is based the behavior of the transformation DJ (Z,, Z2), which is stated below in Theorem 2.2.
1/2 for d N, and Theorem 2.2 (a) Assume that h = k E {1.2}. Define Sa (E 1/2 ({3)) 1/2 assume that 1 = lima (Sa) erists. Then, Tl is finile and
D (1,, Z2) + L' as d + oc. oc.
(b) Assume that + (1,2). Then, Lhk lima D Zn p m
Define Shi == -//2>"1 + > (2*)) 1/2 for d E N, and assume that Lish = lima (trace(S exists. d Then, Lik is finite and
2 (Z, Zz) f. Li as d + 00.
Remark 2.2.1 The structure of the covariance Σ stated in equation (5) imposes some restrictions on the associated constants as stated in part (b) of Theorem 2.2. In particular, the fact that L and T are finite implies that assumption (2) in Theorem 2.1 always holds for the sequence of matrices { } JaEN and (5a JaEN with h, k E (1.2].
Remark 2.2.2 It follows from part (b) of Theorem 2.2 that the statistic we propose is useless for cluster analysis in the homoscedastic case (independently of the difference between μ² and م²²1) because if 2 2 then 1 L'3 Liz. A possibility transformation Dd (12,, 22) increases with Our proposal is use
modify the statistic D²d (z), 22) so that the value of the
D (u, v). :- (12)) v)a|
Σ
(u) with E N. X
Discussion of this transformation, and some numerical results are included in Section I. of Appendix IT.
To simplify notation and avoid technicalities with empty classes, we additionally assume that the observations whose indices belong to the sets C1 f1....,NN N1} and C2 {N1 + 1, I.....N} with N = Nı + N₂ and N1, N2 > 0, were generated by P1 and P2, respectively. In practice, these sets are unknown and in fact our aim is their estimation. We begin with this simplifying assumption for case of notation, and to obtain a clearer exposition of the proposed methodology.
Define the N N matrix I'a, whose (i,j)-th element
N Σ
[D (Z, (Z., Z.) - D (Ze 2,) }
Ta Zz) = 7/3 =
(9)
N
1=1, (.
for S i, S N. Theorem 2.2 and the fact that 1 + j.j in (9) give us the following:
0,
if E Cp. for h 1,2,
d Tag
4
(10)
Thki
if E Cn and Cc, with h # ke E (1,2).
5

---

Published in Transactions on Machine Learning Research (09/2023)
d + Here,
NA
No
N 2 (LS L
N - (L' 12
Note that N1 and N2 are fixed here. Combining the fact stated above in (10), as d -> 20, we obtain
ON, 0
12lNN 17,
La 4 -r==
(11)
721 12iIN,11
0
Let 8 and 3 (for < N) denote the eigenvalues corresponding to the matrices Da and T, respectively. Define the following quantities
N Ka Σ
N
103511 > aa) and Ko = Σ
10 > 0),
(12)
with fastaeN decreasing to 0 as d x (at appropriate rate), and / is the indicator function. The constant Ko clearly equals 2 for the limiting N x N matrix I stated in (11), and hence, correctly identifies the true number of clusters.
Proposition 2.3 Assume N2, 2 1 are fired. Under the assumplions of Theorem 2.2, if L + L's and 13 + Lg then there crists sequence {adldeN CR such that da + 0 and Ka 22 as
This now implies that we can correctly identify the true number clusters asymptotically as d + co. The structure of the matrix I' (11) is straight forward because of the simplifying assumption on the sets C1 and Cz. However, this not a requirement and we will drop it. Proposition 2.3 holds more generally for any permulation of the data points Zi..... Zx- In fact, if the sets Ci and C2 unknown, then the rows/columns of the I matrix will be permuted accordingly. But, the underlying structure remains the same and Proposition 2.3 continues to hold. As a followup of our previous result, any standard clustering method is used on the I' matrix, then we can perfectly cluster all the observations asymptotically (as + oo) because of the structure of the I matrix stated in (11).
prove that if
Definition 2.4 A clustering method can be defined as a map from H to the set {I,...,J) J). Consider the sequence of maps (ValdeN and a second map d. A measure of distance between two clusterings based on the Rand inder (see 847 of Rand (1971)) is defined as follows:
RAN =
1 (2)
Σ ISIGiEN
[[[Wa(z.)) = -a((z)) lb((2; = o(2,)] = =1]
for a fized N 2 2.
Let & be the map which gives the true labels, i.c., o(x,) for j E C₁ with h E (1,2). We can construct a data based vaa by directly applying any clustering technique on the rows or columns of the matrix Ta. Here, we use the k-means algorithm on the rows of Ta-
Mathematically, the k-means algorithm finds J groups (say, ..) with centers CI,...,C. C2 such that o(I1.... -x= E(ix,eG.) |||;- , (2 is minimized. The asymptotic properties of the matrix Ia as d + do (stated above in (11)) imply that differences in the limiting constants should yield perfect clustering. Our next result proves label consistency for this k-means algorithm when 2
Thcorem 2.5 Assume J and m2 O. Further, assume that the conditions in Theorem 2.2 and Proposition 2.3 hold. Then, the clusters will be perfectly identifiable, i.e., BaN 4. 0 as d + xc.
Remark 2.5.1 The well-known Rand index (a measure of similarity) is usually defined as RAN. As a consequence, Theorem 2.5 implies that the usual Rand index goes to d -

---

Published in Transactions on Machine Learning Research (09/2023)
Remark 2.5.2 The structure of the N x N symmetric matrix I stated in (11) continues to hold, and will lead us to perfect clustering for every value of J 2 2 if enough distinct Vi 's exist. In partieular, this happens if we have some h with 1 S h S J such that "Akk * Yhke for every 1 S k + W S J, but other possibilities cxist as well.
Moreover, the procedure described in Proposition 2.3 also works fine with the limit equal to Ko (which the rank T)). One may be tempted to think that it generally coincides with J. But, this is true only for J s 3 and may be different for J 2 4 (as shown in Lemma F.1 of Appendix F). The proof of Lemma F.1 further shows that the condition under which rank(T)) < is quite restrictive. Thus, for simplicity, proposal is to estimate the number of clusters J using Ka (defined in equation (12) above).
2.1.1 Example with GPs
If assume that Σ22 a (for > 0), then we have the following expressions for the scale constants stated in Theorem 2.2:
2a
1+a and L
33
T1 + T2a
+ T2a
+ 032
Thus, it is possible to identify perfectly the clusters as long + 1, since this implies that m22 and 721 both are positive quantities.
2.1.2 Uniform Convergence
In Theorem 2.2, we have proved consistency for finite sets of data points corresponding the transformation D (Z), Z2) defined in (1). We now prove the uniform (on the random sample) convergence of this function as N co. This result will be useful in establishing second result on uniform convergence, which we state in the next section.
Theorem 2.6 Assume the conditions in Theorem 2.2, and let {dr} C N be such that dN - 0o N + co.
Then, we have the following.
a) For h (1,2), let adw (a ,..., ANYT be the eigenvalues of Shxx with dN N. If
log N = o ( 1110
dN
(13)
then it happens that
z
sup 2, Z.EZ
D Dd (Z1. Zz) L1 N N
(14)
b) For any k ke (1,2), let dax = dN
dNyT the eigenvalues of Shkk with dN E N. If
log N =0
dN
(Iaadyllo)
(15)
then it happens that
sup Z, E,EON 2.ECN
r.z Dax (Z1, [2)
Lbk
00 +
(16)
Remark 2.6.1 Assumption (2) holds here, 80
|aayy Ilxe dy dx isisdy
1
dv + 0 as N + x. Thus, if we max 0,
take dy growing fast enough, then it assured that assumptions (13) and (15) hold. The structure of the matrices Sh and Sh for d N with h * le E {1,2] implics that sufficient condition is log o(dy) (also see Proposition H.I in Appendix H).

---

Published in Transactions on Machine Learning Research (09/2023)
3 Transformations with Estimated Gaussian Distributions
In this section, we will discuss the steps to implement the procedure described in Section 2. In practice, the involved distributions and all the associated quantities need to be estimated from the data. Hcre, Z will denote a random element with distribution the GP mixture TT T
Let (t) with & [0, 1] and Xÿ for denote the eigenfunctions and eigenvalues of 23, respectively, We will now make the following assumptions:
A.1 s E(Z <c..
A.2 Assume that 27 > X > ... > satisfying Σj=1 27 < 20.
It is well-known that assumption A.2 implics {o)) forms orthonormal basis of H.
To estimate Σ and its eigenvalues and eigenfunctions, we will use the corresponding empirical quantities. Suppose that we have simple random sample ,....., Zn taken from Pz. Given 8, [0, 1J, we define
N N (.( Zz) (1))
(s,
where (1) = N/ E (t). Consider the corresponding families 17 2 X3 2 ... and Z7 ,_.... of its eigenvalues and eigenvectors, respectively. Note that Σン well as all the 27's and o3 depend on N.
With finite sample, we cannot estimate the entire (infinite) collection of eigenvalues and cigenvectors. Thus, we follow the work of Delaigle and Hall (2012) and Hall and Hosseini-Nasal (2006), and select a non-random decreasing sequence yn going to zero slowly enough as to satisfy limy / n/ = co. We define
R = ij X5 j+1 < ENN 1.
(17)
This definition implies that X3 2 UN for every S R Moreover, we also need that the theoretical eigenvalues are reasonably well separated. To obtain this, given > 0, we further define
R = inf X3 23+1 (1 + 8) - 1.
(18)
We now state empirical analogues of the results stated in Section 2.1.
3.1 Consistency of Clustering
Let Zi...., ZN be simple random sample taken from the GP mixture distribution Pz. Now, Pz and the sets C₁ and C2 (containing information the class labels defined in (6)) are wnknown. Extensions of Thcorems 2.2 and 2.6 to Theorems 3.1 and 3.2 are presented below. The following results will be based on the analysis of the map DR (u, v), which is the transformation D (u, v) defined in (1) with = Åx (defined in (17)), and the pooled covariance matrix 2Z
(sample covariance of which estimated by SZ "RN
RN
the full sample). Here, Day is an abridged notation of D SnN The first result is related to the consistency of the transformation on finite sets.
Ax
Theorem 3.1 Let assumptions A.1 and A.2 and those in Theorem 2.2 hold.
(a) If (1,2), then
Day [Z,
P5 L3 N 00.
(19)
(b) 1f 4 (1,2), then
D Rx Z2) Pg Lis N co.
(20)
8

---

Published in Transactions on Machine Learning Research (09/2023)
We need an increasing sample size in order estimate the parameters consistently. Thus, it is desirable to be able to cluster the increasing number of data points, asymptotically without error. The only way to achieve this is to get some kind of uniform convergence in (19) and (20) when the sample size increases. This is the purpose of Theorem 3.2, which gives us clear evidence that using this transformation would lead to asymptolic perfect separation in the empirical case as well.
Theorem 3.2 Let us assume all the conditions in Theorem 2.6 with log N o((ik) in (18).
(a) For E (1,2), we have that
DAN (Z,, Zz) L/3] 30 as N 00.
sup 21.
(b) For any h E {1,2] with + k, we have that
DAN (Z,, Z2)
L.gh as N
sup ZiEC ZZECN
Remark 3.2.1 Clearly, Theorem 3.1 follows from Theorem 3.2. But, the conditions required for proving the former are weaker and hence, we state it as separate result.
Remark 3.2.2 (Asymptotic perfect identification of clusters) Recall the matrix L' from (9) with d E N. Now, consider the matrix f Aw: which is obtained by replacing sh the matrix I Fx with their estimated values 1B i.c., 16 Fix 6 > O. Theorem 3.2 implies that with probability converging to one as N - 0o, we have
Rx Day (Z,, Zy) with S i $ & N. Define 012 = #1 LL _ + m 1L3 L3
if Z, E Cn for h E [1,2], then (/ s 4e?,
Zg Ch. Zy E Cs for h # E {1,2}, then 1:
V12 S He,
for some H1 > 0. Consequently, if v12 > O, then the elements in I'd will be clustered into well-separated clusters: one around 0 and another one around 812 with probability converging to one.
Similarly, let Pz be a mixture of J(> 2) components and denote
UhA = Th L3 Lake Th L'
L
with 1 S h + h < J. For positive and distinct uhas, the elements in the matrix f'a will be perfectly clustered into + (2) well scparated clustors: one of them around the point and the remaining around the values Uhk (for h < k) with probability converging to one N oc. Therefore, asymplotically, the sequence of matrices { ft. NeN will contain enough information to perfectly cluster all the data points in sample.
3.2 Implementation Issues
We given sample of data points without the labels. Here, consider the N x N estimated matrix fx with the (i,j)-th element as D RN (Z, Z,) (which is just the empirical version of D (Z,, Zj) based on the pooled sample covariance) for 1 < i,j < N and apply any clustering procedure on its rows (or, columms). Note again that we do not need to estimate the unknown constants L/ and L for h, k e {1,2} (stated in Theorem 3.1) for the implementation of our clustoring procedure. The expressions related with K2 and RN are not used in our implementation (see Table 3.2 below for the main differences between theory and practice). In Section 4.2, we give complete details of the implementation of our procedure.

---

Published in Transactions on Machine Learning Research (09/2023)
Table 3.1: Some key differences in our theoretical assumptions and implementations.
Quantity
Theoretical Assumption
Implementation
RN
involves Xys and goes zero as N + 30
estimated using valssaanii clustering (see Wang (2010) for more details)
Ka
involves B and goes to zero as d > 0c
involves 3,8 and estimated using the optishrink function from the R. pack- age denciseR
4 Analysis of Simulated Datasets
For our simulation study, we consider two class problems (J 2). We generated data on a discrete grid of 100 equi-spaced points in the unit interval [0,1] 1] from four different simulation models, which described below. Fix >>6 0.
Define (t) = = (A1/2 Znz+ (un))o()) with E [0, 1] and h = 1,2. Here, the Znys are independent
standard normal (i..., (0, 1)) random variables, (y(() v sin(tjt with t E [0, 1] and 1 s j S 40, Also, Maj = 0 for j > 6, and
set the other components equal to -0.5 1, -0.5, 1, -0.5)* and (0, -0.50.75, 0.75, -0.15, 1.4, 0.1)* for k = 1,2, and \1 = 1/j² and X25 = s/j² for 1 I 40.
This model is from the paper Delaigle and Hall (2012)
II. In this example, X1 ~ B and ~ with pll) = Gl for / E [0, 1] and G G N 4) independent of B. Here, B is the standard Brownian bridge, i.e., a centered GP with Øıj = min(t; (2) t.ty for 4,46 E [0, and i.j E N. Since E[X2(()] E|G fact, the inclusion of μ modifies the covariances because if < t; < t, < 1, then the independence between G and B yields the following:
for í E [0, 1]. the differences in mean never appear in this setting, In
[[X2(t)xx(()]] 4h,t + tz).
This model is from the paper Berrendero et al (2018).
III. Let Xa = ph + E Eny o for h 1,2. Here, Enys are iid. N 1), μ1 and p2(t) with t E [0. 1], X1 j/3 and X2j = /se 1/3 for 1 s 50, and Φx-1 v2 sin(2it)) and Φ2. = 2 cos (2int) for < i < 25 with E [0, 1]. This model is from the paper Dai et al (2017).
IV. This problem consists of two Brownian motions defined in the closed interval [o. 1] with mean fune- tions 1 (t) = 20 (1 () and (t) = 20t(1 - -))--, respectively, for / [o, 1]. For the first class, the eigenfunctions are (t) = 2 (0.5 and associated eigenvalues are \1 = 1/(=(j-0.5)) for 1 و < 15. The second class is similar to the first one, but the eigenvalues are multiplied by / (i.c., X2j / (3/(m(j 0.5)) for S < 15.
This model is from the paper Galeano al (2015).
We set s = 1 for location only problems. In location and scale problems, we fixed s = 3, while for scale only problems the mean functions z and 2 were set be the constant function and = 3 was retained.
4.1 Choice of d
A critical issuc here is selection of the optimal dimension of the projected space for a given set of data points (i.e., a fixed value of (N1, N2) or N). Let us recall Theorem 3.1. According this result, we expect the rows of the matrix Pdd to form two clearly separated clusters depending the class label of the observation for
10

---

Published in Transactions on Machine Learning Research (09/2023)
large values of d. To demonstrate this, construct sequence of images which shows how this separation varies with increasing values of d.
We generated samples of size 250 from each of the two classes for the 'scale case' of Example II. For the purpose demonstration, the first 250 observations correspond to the first GD, while the next 250 observations the second. Figure below shows the heatmap for increasing values of d, and we observe the best concentration d 80. However, some noise in the off-diagonal submatrices for d 80 (compared to d 60) makes to consider that the optimum could be somewhere between the values 60 and 80.
(a) d= 20
(b) d=
(c) 60
(d) d 80
Figure/ 1: Heatmap of f'a for varying values of d.
Clearly, the choice d quite important as is the dimension the subspace where we project our obsorvations (for fixed sample size). We observe from Figure that its estimation quite crucial. The next subsection contains further details on the choice of d.
4.2 Clustering Procedure Practice
To implement the clustering method, we usee vros-vvaliattin (CV) choose the dimension d suitably. We use the idea developed by Wang (2010), which state briefly here. Given B E N, split the data into three random subsets (say, S16, S2t and Sys) cach of cqual size for 1 < b S B. For each value of b, treat the data points in SIa and S25 the training sets, and S3, the validation set. For a fixed value of d and given a clustering algorithm, the two training sets S1 and S25 are used to construct two cluster assignments. An appropriate distance between these two cluster assigments (say, D) is computed based on the validation set Syt (see Section 2 of Wang (2010) for more details). We repeat this partitioning BB= 50) times and average it over these B samples to get LV Define dev = arg minggaeN DCV. A pseudo-code for this procedure is given Algorithm 2 (sec Section O of the Appendix).
Recall the structure of the I matrix stated (11), and also see Figure As mentioned Section 3.2, the number of clusters were estimated using the quantity Ka described in Section 2.1 (see (12)). To implement the procedure in practice, one needs to estimate the sequence {au]deN- We have used the function optishrink available in the R package denoiseR. This function extracts a low-rank signal from Gaussian noisy dala using the optimal shrinker of singular values. The low rank structure of the I' matrix motivates us to directly apply this function on Pa. A pseudo-code (alongwith stepwise computational complexity) of our methods is given in Algorithm 1 (see Section
the Appendix).
Our overall implementation yielded quite desirable results in our numorical study (sce Tables 4.1 and 4.2 below). We can apply any clustering method the transformed data Пa. In addition to the k-means algorithm (CD-k-means) discussed in Theorem 2.5, considered spectral clustering (CD+spectral) and Gaussian mixture models (CD-mclust).. Here, the acronym CD corresponds to our proposed transformation. One may refer to the book by Hastie et al (2009) for details on these three popular clustering methods. The R codes for our methods are available here: GP clustering.
We considered several methods for comparison. The first three methods are those that we have used on the matrix fa but applied directly to functional data, namely, the k-means algorithm (k-means), spectral clustering (spectral) and Gaussian mixture models (melust). Several competent methods for functional
11

---

Published in Transactions on Machine Learning Research (09/2023)
clustering using functional mixed mixture models are implemented in the function funcit from the R package funcy. We report this method as funclust. The methodology developed by Chiou and Li (2007) is available in the function FClust from the R package fdaspace using two clustering techniques 'EMcluster' (CL1) and "CFCC" (CL2). We have reported the minimum value, and stated it as CL. In Delaigle et al (2019), the authors developed functional clustering based on the k-means using basis functions. We implemented this method for two choices of the basis functions, namely, Haar and PC, and reported the best result among these two (we call it DHP). We have not used the DB2 basis for our comparisons because it requires the grid points to be of power of 2. The DHP method is available from the journal website, and used those Matlab codes for our comparisons.
We conducted simulations based on models to IV, which were introduced in the beginming of Section 4. We did not consider the location only scenario as our proposed method useless in such cases (recall part (c) of Theorem 2.2). However, we have some discussion and additional results in Section L in Appendix II for this scenario. The sample size of each class was set to be 250. Our experiment was replicated 100 times, and the results are reported in Tables 4.1 and 4.2 below. To measure the similarity between two cluster assignments, computed the adjusted Rand index using the function RRand in the R package phyclust. One minus the adjusted Rand index (we call it adjusted Rand distance) is reported in tables below, where the minimum is marked in bold and the second lowest is in italics.
It is worth noting that all the competing methods require the number clusters as an input variable, and we have run these methods with
(the true number of clusters). However, when applying the CD procedure we have estimated the number of clusters following the procedure described above. We obtained the correct value in more than 99% of the cases (across all four examples for both scenarios) in our simulation study.
Table 4.1: Adjusted Rand distances for different GPs with differences in locations and scales (with standard error in brackets).
CD
Ex.
k-means
spectral
mclust
funclust
CL
DHP
k-means
spectral
mclust
0.0632
0.0280
1.0000
0.1541
0.0239
0.0818
0.0654
0. 0.0028
0.0001
(0.0007)
(0.0016)
(0.0000)
(0.0017)
(0.0007)
(0.0025)
(0.0010)
(0.0002)
(0.0001)
II
0.5377
0.4210
0.9700
0.8222
0.5767
0.5149
0.2104
0.3464
0.0738
(0.0046)
(0.0042)
(0.0017)
(0.0027)
(0.0045)
(0.0049)
(0.0030)
(0.0045)
(0.0014)
III
0.4250
0.4460
1.0000
0.3858
0.2891
0.4137
0.1367
0.2596
0.0250
(0.0017)
(0.0048)
(0.0000)
(0.0003)
(0.0000)
(0.0054)
(0.0024)
(0.0042)
(0.0016)
TV
0.4945
0.3788
1.0000
0.3975
0.1833
0.1379
0.0316
0.0000
0.0000
(0.0056)
(0.0047)
(0.0000)
(0.0011)
(0.0000)
(0.0033)
(0.0001)
(0.0000)
(0.0000)
(a) Fx.
(b) Fx.
(c) Ex.
(d) Fx.
Figure 2: Representative curves for the four examples having differences in locations and scales. The mean curves marked in bold, while the two classes marked in red and blue colors, respectively. Please zoom in for better resolution.
In the first setting, we considered clustering problems with differences their location and scale param- eters. Class specific representative plots are given Figure 2. Usefulness of the proposed transformation is clear from Table 4.1. Our method attains the first position across all examples, while in Example TV
12

---

Published in Transactions on Machine Learning Research (09/2023)
obtain perfect clustering. Although there is location difference in Example II, sub-optimal performance of our method is probably due to low signal from the difference between the two covariance structures. CT. attains the second best performance in the first three examples among the competing methods, while DHP performs better than CL in Example IV.
In the next setting, we dealt with differences only in scale parameters. Class specific representative plots are given in Figure 3. It is clear from Table 4.2 that the separation in scatters is captured very well by the proposed transformation F'. Moreover, improvement in Example II compared to Table 4.1). In Example II, all (except possibly, funclust) the competing methods perform quite poorly across all examples. However, its performance is far below our method CD even in this case. The performances of k-means and DHP are similar, and quite bad in this scenario. Generally, the results Table 4.2 suggest that existing methods fail to judiciously capture information if present only in the scale parameters.
method again leads to perfect clustering (with significant
Table 4.2: Adjusted Rand distances for different GPs with differences only in scales (with standard error in brackets).
CD
Ex.
k-means
spectral
mclust
funclust
CL
DHP
k-means
spectral
mclust
0.9900
0.9992
1.0000
0.9776
0.8269
0.9966
0.0063
0.0001
0.0000
(0.0014)
(0.0007)
(0.0000)
(0.0006)
(0.0000)
(0.0005)
(0.0002)
(0.0001)
(0.0000)
11
0.9985
0.9962
0.9986
0.5004
0.9065
0.9999
0.0091
0.0100
0.0084
(0,0005)
(0.0007)
(0.0005)
(0.0049)
(0.0007)
(0.0003)
(0.0036)
(0.0047)
(0.0019)
IIT
0.9944
0.9990
1.0000
0.9956
0.9994
0.9967
0.1789
0.2872
0.0352
(0.0009)
(0.0001)
(0.0000)
(0.0000)
(0.0000)
(0.0007)
(0.0026)
(0.0032)
(0.0007)
IV
0.9927
0.9984
1.0000
1.0006
0.8464
0.9980
0.0102
0.0014
0.0005
(0.0002)
(0.0005)
(0.0000)
(0.0000)
(0.0000)
(0.0006)
(0.0021)
(0.0013)
(0.0004)
(a) Ex.
(b) Exc. 2
(c) Ex.
(d) Ex. 4
Figure Representative curves for the four examples having differences only in scales. The mean curves marked in bold, while the two classes marked in red and blue colors, respectively. Please zoom in for better resolution.
After applying the transformation f'a, we had used three methods for clustering the transformed observations (see Tables 4.1 and 4.2). Additionally, when these methods are applied directly to the functional data, we observe that k-means and spectral perform fairly well, but their performances deteriorate sharply for the case with differences only in the scales. The mclust algorithm performs worst (possibly due to the presence of low signal in the locations) among the three clustering methods, while CD-mclust clearly achieves the best results with substantial improvement (see Table 4.1). It is worth mentioning that the Rand distances improve substantially for all three methods when applied on the transformed matrix f'd
All the three usual clustering methods seem useless (possibly due to the presence of signal only in the scales) from the results in Table 4.2, while their CD counterparts clearly lead to excellent performances with CD- mclust yielding the overall best result. Overall, the best performance of this method may be attributed the fact that the signal gets amplified the transformed matrix and the specific structure (recall equation (11)) of this matrix.
13

---

Published in Transactions on Machine Learning Research (09/2023)
5 Analysis of Benchmark Datasets
We have applied our proposed methods to some benchmark data scts, Wheat (from the R pack- age fds), Satellite (available at https://www math. univ-toulouse. f /index.. html), Cars (kindly provided by the first author of Torrecilla et al (2020)) and Velib (from the R package funFEM).
To evaluate the clustering algorithms, we ran a single execution (without splitting). Class assigments are already available for the Wheat dataset. The Satellite data has been analyzed in detail in the paper Dabo- Niang et al (2007), where the authors split the curves into two clusters 'unimodal' and 'multimodal,, The authors of this paper kindly shared the exact cluster assignments for this data set with us. The Cars data contains asset log-returns of the car companies Tesla, General Motors and BMW (see Torrecilla et al (2020) for more details). However, the rank of the estimated f matrix was two for this data set, and our mothod detected only two distinct clusters. This coherent with Torrecilla et (2020), where the authors had noted that assets of General Motors and BMW were very similar and quite difficult to distinguish. So, we merged General Motors with BMW while assigning the class labels for this data set. Consequently, the number of clusters was set to be two for all the competing methods. The Velib data analyzed by Bouveyron et al (2015), where the authors identified the optimal number of clusters to be ten using the funFEM algorithm. Setting J = 10 (this information was provided to the competing methods), we determined the class labels of the observations using this algorithm. We report the adjusted Rand distance for these four data sets in Table 5.1. Competitive performance of our proposed methodology w.r.t. the competing methods is clear from the results given below.
Table 5.1: Adjusted Rand distances for different clustering methods.
CD
Data
J
N
d
k-means
spectral
mclust
funclust
CL
DHP
k-means
spectral
mclust
Wheat
2
100
701
0.6960
0.7860
0.6978
0,6960
0.8058
0.5730
0.7859
0.8057
0.3644
Satellite
472
70
0.6072
0.8947
0.3525
0.6072
0.6060
0.7253
0.8861
0.9954
0.4448
Cars
2
90
32
0.8856
0.9230
0.8834
0.8856
0.9650
0.9088
0.5385
0.9313
0.4680
Velib
10
1189
181
0.3755
0.5464
0.8561
0.6738
0.5672
0.5872
0.5408
* Rt package is now archived; ** Matlab code available only for two classes.
To get a better understanding of the performance of our proposed method, we further computed the well- known average purity function correspoinding to the CD) method having the minimum adjusted Rand dis- tance. A value of average purity function close to one indicates good performance of a method. We obtained the values as 0.9000, 0.8622, 0.8666 and 0.6038 for the Wheat data, the Satellite data, the Cars data and the Velib data, respectively. Overall, our proposed method CD yields quite promising results in all four benchmark data sets, with CD-melust yielding the most stable performance.
Acknowledgements
The first author has been partially supported by grants PID2021-128314ABB-10 and PID2022-1332337BBBBB UE and "ERDF A way of making Europe'. Both
funded by MCIN/AEI, 1 the authors would like to thank the Action Editor Prof. Brian Kulis and three anonymous reviewers for their constructive comments and suggestions that improved the paper.
14

---

Published in Transactions on Machine Learning Research (09/2023)
References
Baillo, A., Cuevas, A. and Cuesta-Albertos, J. A. (2011) Supervised classification for family of Gaussian functional models. Scand. J. Statist., Vol. 38, 480-498.
Berrendero, J. R., Cuevas, A. and Torrecilla, J. L. (2018) On the use of reproducing kernel Hilbert spaces in functional classification. J. Amer. Stat. Assoc., Vol. 113, 1210-1218.
Bouveyron, C., Cóme, E. and Jacques, J. (2015) The discriminative functional mixture model for compar- ative analysis of bike sharing systems. Ann. Appl. Statist., Vol. 9, pp. 1726-1760.
Chiou, J. M., and Li, P. L. (2007) Functional clustering and identifying substructures of longitudinal data. J. R. Statist. Soc. B, Vol. 69, 679-699.
Cuesta-Albertos, J.A. and Dutta, S. (2022) On perfect classification and clustering for Gaussian processes. arXiv preprint,
Dabo-Niang, S., Ferraty, F., and Vieu, P. (2007) On the using of modal curves for radar waveforms classifi- cation. Comput. Statist. Data Anal, Vol. 51(10), 4878-4890.
Dai, X., Muller, H-G. and Yao, F. (2017) Optimal Bayes classifiers for functional data and density ratios. Biometrika, Vol. 104, 545-560.
Delaigle, A. and Hall, P. (2012) Achieving near perfect classification for functional data. J. R. Stalist, Soc.
B, Vol. 74, 267-286.
Delaigle, A. and Hall, P. (2013) Classification using censored functional data. J. Amer. Stat. Assoc., Vol. 108, 1269-1283.
Delaigle, A., Hall, P. and Pham, T. (2019) Clustering functional data into groups using projections. J. R. Statist. Soc. B, Vol. 81, 271-304.
Feldman, J. (1958) Equivalence and perpendicularity of Gaussian processes. Pacific J. Math., Vol. 8, No. 4, 699-708.
Galeano, P., Joseph, E. and Lillo, R. E. (2015) The Mahalanobis distance for functional data with applications to classification. Technometrics, Vol. 57, 281-291.
Hajek, J. (1958). A property of J-divergences of marginal probability distributions. Cz. Math. J., Vol. 8, 460-462.
Hall, P. and Hosseini-Nases M. (2006) On properties of functional principal components analysis. J. R. Statist. Soc. B, Vol. 68, 109-126.
Hastie, T., Tibshirani, R. and Friedman, J. (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. New York: Springer.
Laurent, B. and Massart, P. (2000) Adaptive estimation of quadratic functional by model selection. Ann. Statist., 28, 1302 1338
Rand, W. (1971) Objective criteria for the evaluation of clustering methods. J. Amer. Stot. Assoc., 66:336, 846-850.
Rao, C. R. and Varadarajan, V. S. (1963) Discrimination of Gaussian processes. Sankhya, Series A, Vol. 25, 303-330.
Shepp, L. A. (1966b) Gaussian measures in function space. Pacific J. Math., Vol. 17, 167-173.
Torrecilla, J. L., Ramos- Carreño, C., Sänchez-Monttaiie M. and Suárez, A. (2020) Optimal classification of Gaussian processes in homo- and heteroscedastic settings. Statist. Comput., Vol. 30, 1091-1111.
Wang, (2010) Consistent selection of the number of clusters via cross-validation. Biometrika, Vol. 97, 893-904.
15

---

Published in Transactions on Machine Learning Research (09/2023)
Appendix I: Proofs and Mathematical Details
A Proof of Theorem 2.1
Fix d E N. The d-dimensional random vector (Z b)a has Gaussian distribution with mean equal to (n b)a and covariance matrix equal to Σg- Now, [I(As) b)all is equal to the square of the norm of a d-dimensional normal variable with mean m. = (Aa)-1/2(p b)a and covariance matrix Sa = Therefore, if ud is a d-dimensional vector with centered normal distribution and covariance matrix equal to Sa, then
(Aa)
D (z. ~ 1 (ma + un ma + ua)
1 d (|mmal + ||uall + 2(ma. ua))
(21)
By assumption (3), we have
lim d d ||maal Tap.
Let consider the second term in (21). Fix basis in Va spanned by the eigenvectors of Sa- Note that this term is not dependent on Ly- Denote uı (ud,1s.... ua and ma (m......... ma,a) in this basis. Therefore, the random variables (ua,i)² with 1 S i < d are independent with means equal to of for 1 i S d and Ei= (ua,a)² ~ Eiz (u,).. Here, {ushsssa is sequence of independent and identically distributed (i.i.d.) real variables with the standard normal distribution. We split the proof into two cases.
A.0.1 Lg is finite
Fix c > 0. Taking into account that the variance of x² distribution with one degree of freedom two and using Chebychev's 's inequality, we have that
P
|ual trace(Sa)) Ze
= P
((ua))
aaf
26
d
1=1
2
2
2
11
0
2
which converges to zero by assumptions (4) and (2). Consequently, we have shown that
1 uu 1| ttrace(Sa) f 0 as d oc,
and assumption (4) gives
d ||al| E, L.s as d + oc.
A.0.2 Lg is infinite
We have that
d
P
((ua,)) - a.)
26 = P
(lu:) -1) 1)
2
ad
i=
i=1
i=
d
2
a
d
i=1
2 ||
2 Σ
1=

---

Published in Transactions on Machine Learning Research (09/2023)
which converges to zero because Ls oo and assumption (2). Thus, have shown that
2 oj (d [ traco(Sa)) f 0 as d oc.
(22)
Consequently, 2 /|ll converges to co at the same rate as ttraco(Sa).
Concerning the last term (21), we have (ma, ua) E ma, We split the proof into cases.
A.0.3 by is finite
Fix e > O, and define a² = (a))) Using Chebychev's inequality again, we get
mm uaحl > e S ddd (max)) S 2dd la |mall"
which converges to zero by assumptions (3) and (2), and the proposition is proved in this case.
A.0.4 Iu infinite
The result follows from equation (21) and the previous results, if we are able to show that the sequence of real-valued random variables
(m us)
Wd max mall 13, |ua]
converges to zero in probability as d - co. In turn, this will be fixed if we show that every subsequence of {wa} contains a new su subsequence which satisfies this property. Thus, let [wa,} be a subsequence of {wa)
and let us consider the associated subsequences {[|ma |l} and [|] ]] Obviously, there exists a further subsequence {dg-) such that one of the following holds:
||ma,- 112
(i) limd, Srace(Sdlz.
|| 00.
(ii) limde trace(Sdr.)
|I 112 C.
(iii) There exists finite C > 0 such that litma tracel (Sc.)
Note that in cases (i) and (iii), we have Ls co. To simplify notation, we denote the sequence {Sd,1 by {S,), and similarly for the remaining ones. In case (i), since equation (22) shows that
| Pg as + oc. (Ta
(23)
mm/l |
have
P as h - co. Consequently,
(m,, un)? |
Imal
lim [wyl
lim h
< lim
in probability.
h
h
|u
If (ii) holds, have that |wal <
u. my
Since
E || || trace(Sa), we have that
U 1
and also, in P.
my
this case wh 4 0 as + oo.
case (iii), taking into account that equation (23) now holds, it is enough to show that
(mg, un) P 0 as h + oc,
(Strcee(p)
2

---

Published in Transactions on Machine Learning Research (09/2023)
Fix > O. We have that
(m, un) C'raaeeSS
h
2
||a loc
mm 1|
>(
3 .
c
a
C
2i=1 a
trace(Sn)
which converges to zero by assumptions (4) and (2).
B On Assumptions (4) and (2)
The next lemma shows that if Ls co, then assumption (4) implics assumption (2).
Lemma B.1 Let faata21 be sequence of real positive numbers such that lim à Ei= a, erists, and is finite. Then, it happens that lima 1loa 0.
Proof: Fix d E N, and denote Aa = Σi= dt. We have that
ad
Ad
Ad-1 d d-1 d-1 d
d d
and consequently, 0 == lima 3 Given >> 0, there exists di > 0 such that if d d1, then o S and d2 2 di such that
a
sup < c. IsiSd d2
Let > d2 and take 1 i S d. So, we have that if < d1, then x < dz This completes the proof.
se and if > d1, then 4 s Se.
C Proof of Theorem 2.2
First, note that (Z)- Z2)a is d-dimensional normal vector, with mean (μ 12 and covariance 2 7 for h,k E (1.2}. Now, the proof will based in the following lemma.
Lemma C.1 Assume that 1s (1,2)- Define Sh St and Shi in Theorem 2.2. Then,
i) L. = lima D 27 (n ==)= =0, and
i) lim sups Atrace(SS) < 0o, lim lim (trace(Sh) oc and lim sup, (trace(Shh) oc.
Proof: Let us denote ===mp>>> 22 + T252e, (122 - 2e) and 722 = क1 From (5), we have that
D = Ex + T12H4A1
From here, the Sherman-Morrson formula gives
112 (2)) hand (E
(Ez)-1 -()==(E))-- (Et)-1
1 + mm (Ez)
Since (Ea) is positive definite (p.d.) for all E N, this now implies that.
T22 -E 1+ T22 (n) (ET) pc) + TI2H (Et)
1 (E Hd <
1
0 h Tl (> pd (1
pd
122
and the proof that I h.k 0 trivially ends from definition of T.
3

---

Published in Transactions on Machine Learning Research (09/2023)
To handle the terms L (L's is identical) and L recall the Woodbury matrix identity:
(U - CUV-V11 -U+UV-
Using this identity, we have
- (=
E2 )-1
Ba,
where Bd (TT 2 x/=2 ( + (13Hap3))->>
If and V are p.d., then VU is p.d. In Bd, both the matrices Σd and (T2Id + =12/4a47) are symmetric and p.d., and this implies that Ba is also p.d. Further, (>/ and Ba are p.d. which now implies that (12 (e( )1/2 is p.d. Recall that trace is linear map. Now,
trace(S))
2 ttrace (" 222 (22)
2 2trace (trace(L/TTLa) trace 12)) Ba)
2 ttrace(L/TTla) trace (acee(22) ((> /1/2)///233 Ba (a(2a))11))
trace(2/m1l) 2d/MT1
We have proved that lim supa (trace(S2)) < 2/ Similarly, it is proved that lim supa (rrace(S)) <2/T2. Fi- nally, we have
trace (S Stak)
= trace
-1/25/1/27 - (=2)) 11/2) + trace((E2)) 1/25. - (
d
d
T 1
from where we conclude that lim sup, Itracee Sht) < 1/(
- -122 = O₂ and Σ + >> =" we take Aa = S according Remark 2.1.2, Lemma C.1 gives that assumption (2) holds for this selection of Ag- Therefore, (7) follows from Theorem 2.1 because in this case I and once assume that L3 exists, Lemma C.1 gives that it is finite.
To prove part (a) in Theorem 2.2, we have h = So, (12m-122
In case (b), we have h * k. We take Ad Shat and b = n² z Similarly, as in (a), we have that assumption (2) also holds in this case and Theorem 2.1 implies
D (Z,, Z2) F. Lhk + LS d oo.
Now, (8) follows because Lemma C.1 gives that Luk 0 and, also, that Lish is finite.
Proof of Proposition 2.3
Under the conditions of Proposition 2.3, the number of significant (unique) cigenvalues of the matrix I Recall that N is fixed here.
Consider the standardized distance matrix Da with the (i.j)-th element as D (Z. Zj) for 1 5 i,j N and d E N. We have a sequence of matrices Da & Do d + (componentwise). Since the map D to I' is clearly continuous w.r.t. this convergence, we have that l' P. l' as d + oc. Let us denote the eigenvalues of Pa (respectively, I) to be $ (respectively, $1.... 3x). ). Since cigenvalues are continuous functions of the respective matrices, we have 34 P. 3g as + xo for all < S.
Let us now look into the following:
Σ Ie > ad) 5 N
N
I(159) > 0) as + 00
4

---

Published in Transactions on Machine Learning Research (09/2023)
with Pd + 0 as d - at appropriate rate. Recall that the limiting quantity the right should give us the correct number of clusters. Consider the sequence (1/m)meN- Let us take i E (l...., N] such that 3 4 0 as d + co. Thus, for every e, 5 > 0 there exists Da_ such that if d 2 Da, then
P||Sy| > 81 <.
In particular, if we take å = é 1/m, there exists D', such that if d 2 D then we have
1
|
<
m
Without loss of generality, we can assume that Dj D2 .., and consider the sequence
2
if < i < Dj.
d
1 if D'm Si < Di for some m 2 1.
Then, obviously a² - and
P [T(80T) a%) > 0] = [| > a²]] az.
If we define as sup(a's 0 = 0}, and i satisfies that 3 0, then I(1,,1 > aa) 0 d - oc. A similar
reasoning allows us also to conclude that if |, > 0, then I(18 > ad) L as d + 00.
E Proof of Theorem 2.5
In this proof, we use the superindex d in 9 to emphasize that the groupings can change with the dimension d N. Proposition 2.3 implies that Ka = 2 with probability converging to one.
Note that o( has an alternative mathematical expression as
E
Σ
Σ A=1 u.vE
u v/
(24)
2
where |G] denotes the cardinality of the set G. Let us denote the rows/columns of Pd ,,..., 2x. The structure of rd implies that ll ~ & 0 as d + oo iff i,j E Cn for h E [1,2]. So, if each G for h 1,2 contains observations from the same population, then da 9 P 0 as d
Let assume that the contrary, there exists subsequence of dimensions {dx} such that for every k there exists at least a couple of points is with is E g and ja G (say). Since the number of points is finite, there exists a further subsequence {dg-) } such that both sequences (ix) and (ik- are constant. Therefore, for those subsequences, we will have
lim inf 0a(91, 2 lim My
du 1|2 -0-
722
712 >
So, for the minimization ( 9h),, each 9 must contain all observations from a single population with probability converging to one as the dimension increases. This proves the convergence in probability of the Rand index Ra,N zero as + 2o.
F Rank of the Matrix r
Identifying number of clusters from the matrix T is equivalent to finding the rank of the matrix T.
5

---

Published in Transactions on Machine Learning Research (09/2023)
Lemma F.1 The rank of the matrix I is less then or equal to J. Moreover, equality is guaranteed only when 3.
Proof: Trivially, rank(I)) S J. Let denote the reduced Echelon form of N x N matrix I' as 1 Thus, the matrix O is J x J symmetric matric with 7ij > 0 and distinct when i A J, while Tii 0.
Moreover, for J = 3, we have
712
713
det(TT) det
712
0
723
722 + 0.
V13
Y23
0
In the case J = 4, if 7/2 = Y13724 + 714723 + 2 713371442377 723722
then simple computation gives that det(TT) = O. This happens, for instance, we consider the following matrix (with all positive and dis- linel off-diagonal entries):
734
1
2
0
3
4
3
0
5
4
5
0
where / 4 0.
G Proof of Theorem 2.6
In order to simplify the writing, we will write d instead of dx- We will use the notation ||aa||2 == The real r.v.s {u,} are assumed to be iid. with standard normal distribution.
(E (Ei)( 1/2
The following lemma is deduced from Lemma in Laurent and Massart (2000) on p. 1325, after some simple computations, laking into account that ||aall2 ||aallee We state it here for further reference.
Lemma G.1 If Za Eie (u - 1) and x 2 1, then
Zal 2 4x 4z||oall><) s 2cxp(-xx)
We will also employ the following well known bound for the tail of the standard normal distribution:
2 2 exp(-t8/2 for all t2 1.
PINN 2) 8
(25)
Proof of Theorem 2.6: Let us show part b). The proof of (14) is similar that of (16). We use the notation ma 112) (n 2=22) )a and u (E2))1/2(Z,, 44)a with d N, where Z; is generic observation with distribution P. for 1, 2. Moreover, with obvious abuse of notation, we will often write u' E CN
with d E N for 2.
Recall that Ln = and Ls < oo (see part (c) in Theorem 2.2). Repeating the first steps in the proof of Theorem 2.1, we have that
rz Ded (Z), Zz)
z trace(S1?
sup 2 ECN 2eECN
1
m
sup
1/ u
trace (11)
(26)
u ele EC
+ 2 sup
(m u uz)l
(27)
u ECN eee

---

Published in Transactions on Machine Learning Research (09/2023)
and it is enough to prove that the terms in (26) and (27) converge to zero in probability.
The first term in (26) couverges to zero by first part of (c) in Theorem 2.2. Concerning the second term, let Ni, N2 be the number of elements in CN and CØ respectively. Since Ni + N2 = N, it is clear that N1 & /4/4. Let E 0. We have that
PN = P
sup ef xPECy
d |l uîll trace(()))
>
{ {1211 [|llu - --2ll
traco (512)
>
Lu JPEEN
N²
1 u²il 1---ttan trace((2)))
P
(28)
where u² and 2 are associated with some Z1 E CN and Z2 E C, respectively. However, clear that
2 1 uill - - trace
1
Σ
(4 1).
==1
Take x = c4/(aalal)))) By assumption (15), we have d/ |aaaloo + and eventually x 2 1. So, from Lemma
G.1, we obtain
N PN S
P
Σ of (u 1) ed S exp
1
1|| (Talaad)
(-1 Ed
+2log
log logg))
4
which converges to zero by assumption (15).
Fix >> 0. For the third term, in equation (27) we have that
PN
P
sup
(ma, ut -uzil
EC EPE
2
P
Mmd u
uzحl
4
N
mai mai(a (a/)
4
d
>E
N²
d
P
[N(0, 1)1
>
(ma Pay
2 d
23/2 nT 1/2 exp
+ 22log N Lisa(ma))
(29)
2²
d 2 |1a Eisa(mai) + 2 22
1/22 exp
which converges to 0 because the fact that I, = 0 (see part (c) in Theorem 2.2) and (15). The same assumption allows us to apply inequality (25) equation (29).
H Result Related to Remark 2.6.1
log N
Proposition H.1 Under assumptions of Theorem 2.2, if we assume that and (15) hold.
- 0, then conditions (13)
dx

---

Published in Transactions on Machine Learning Research (09/2023)
Proof: Fix h (1, and note that
(22)-1 (Th Zz " Pa.
where Pa 2d + TA Z. (T2) nT 2- -1 p.d. matrix and T3 the matrix mg.a + TTTTL1E
112 - lddl 7T with 1 + h E (1. Further,
- 22221
La + 2
Th
From here, Weyl's incquality gives
S Ommin
1
2 7t ) Ommin 24 ( (>2
1
(30)
Th
Th
Note the fact that the eigenvalues of the matrices AB and BA are same. So, the matrices S and > (Ez)-1 will have the same eigenvalues. Furthermore, the eigenvalues of Sy are the inverses of the eigenvalues of Σ (EZ")) Thus, (30) gives that
2 Th
Omaz (S < (free of d).
(31)
We now have
dy dw
01 dN log N
log N
0 as N 00.
dy
log N dN
Equation (31) now implics that condition (13) holds we assume
+ as N 00.
Fix h + k (1,2). Our second matrix of interest is
S = 112(-11) + M2 (E3)
Since the matrices are symmetric, have
Ormaz (h S amaz (5
1/2 Zn dd
(EZ)- -1/2) + Omaan (5%))
1/2 1/ Zu 12/)-1//3
Again, the eigenvalues of (Σ7) 1/2
Z.
(E7)
1/2 and of Σz (E will be equal for i h K. So,
aman (5taht <
Oman (2))) (MZ)-1 + Omaz (2) (E2) -1)
amin (> (2 1 aman 22 (2 -
(using equation (30)).
Th
log 0 implies (15) holds.
From here, similarly as before, we would obtain that
dx
Proof of Theorem 3.1
Recall that subspaces Va generated by estimates of the first d eigenfunctions of the covariance of Z.
We begin with some notation and preliminary results which have been taken from Delaigle and Hall (2012) and Hall and Hosseini-Nasee (2006), or follow directly from the results there. Then, we will give the proof of Theorem 3.1. For every n E N, let us consider
A (1) (2)
57 = min(AA AA= \et1).
8

---

Published in Transactions on Machine Learning Research (09/2023)
In Delaigle and Hall (2012) and Hall and Hosscini-Nasale (2006), it is shown that if j 2 1, then
P2j X,l - Âz,
(32)
and that if j S R2 (recall the definition of R in (17)), then
| Q51 v 8 (00))((0))-
(33)
Åz = Op(NV112).
(34)
R roand RN S 1Z
(35)
Moreover, if S R there exists k S j such that
og = > Ak+1 2 A i441 2åz AN 2åz 7n op(Πx),
(36)
where we have applied (32) and (17). Using (34) and the assumption on V|N, we can conclude that V/N from an index onward. Thus, (36) and (33) yield
2Âz
Âz
| - 1
(37)
2Âz
UN
From (32), (17) and (34), we obtain that
A 2 Åz 2 nN
Ô = 1/ + +p(NN).
(38)
Now, we are in position to prove Theorem 3.1.
Proof of Theorem 3.1: The proof is based on Lemma L.1. The result follows trivially from this lemma, the fact that R2 f do as n + oo, and the result in Theorem 2.2.
Lemma I.1 Under the assumptions in Theorem 3.1, happens that
[DR (Zi, Z2) - Day (Z, Zz) f 0 as n + cc.
Proof. For fixed Z1, Zz, let us denote u Z1 Zz. Obviously, ||u|| = O(1) a.s. Let us denote (uI...., " " ) and (01...., 22z) )" to be the projections of u on the subspaces generated by the first A% eigenvectors the respectively, when written in the basis generated by those eigenvectors. Now, define
matrices Σ and Sz 2 12
a -u, 33) J u(t)o)(()))) (t)dt
and similarly define u, for E N. Fix N and take j 5 R2 We now have
(u,)² A2 S
(b,)²
új
Uj
+ úg
=
(2))
(2
(Az)
(2y)1/2
Uy -リー
1/2 (2,)
Uj
+ ú
+ ûg
1//
(1,2,)2/2
(A2)
(2,)112
We analyze each term in this expression separately as follows:
U - 0g
(22) / [u()))l,(())-8)) &
5
(2,)
||u|| 110 - ا11
(22)
Âz
< 8//2||ul
1/2 - 2åz)
< 81/2 ||/2|||| Áz (mm
33/2 + or (1r 33/2)),
(39)
9

---

Published in Transactions on Machine Learning Research (09/2023)
where we have applied the Cauchy-Schwand inequality, (37), (34) and (38). On the other hand, have
0 (4y) (2,))//2
L |u(t)l/;(()dd (2)) + (2,)1/2) (212))1/2
1A -
(Az2,)1/2
Åz
||u|| (2)) + (2,) (A12,)1/2
2z
| 3/2 + (1m3321))
(40)
where we have applied (17) and (38). Concerning the final term, using (38) and (17) again, we obtain that
uj (2,)1/2
+ ū
s ||u|| (Ajjuz + (2,)1/2 S || (1m 1/2 + op(m 1/2).
(41)
(2,)1/2
Now, define C = -8//2 + 1. Combining (39), (40), (41); (35) and (34), we get the following:
(u,)
(a,)²
D ng (2. D2 Z Zz)
R
Aj
j=1
C'||u|| (nm² + opp(m
Op(n' -1/2 0p(mm112mmm
By construction, The is such that any, + co. So, we have 1Dzz (Z1, Σn) - D Daz (Z1, Z2)1 - 0 as - 0o, and the lemma is proved.
J Proof of Theorem 3.2
We now prove the following lemma.
Lemma J.1 Under the assumptions in Theorem 3.2, we have that PpRN 2 RN] + 1 as N 00.
Proof: Fix N E N. From (32), we have that
inf (2, - )5-1) inf 1 RN
(A, A9t1) 2å2 2 (1+8)TN 2Âz.
RS
and the proof ends because (34) and the fact that UN 2 N-1 imply that P(önN 2Ωz 2 0] - 1.
In this setting, recall that L, 0 and Ls oo (see part (b) of Theorem 2.2). We will only prove part b); part a) being similar. W.l.ogg we will assume that h = and k = 2. Recall that for every Z and Zz. we have
DRN (2, Zz) = 1
AN Z
Ze, o,}² and D (Z,, Uz) =
1
RN Z
Zz, S2
RN
Aj
RN
j=1
j=1
We are going consider the function
1
(Z Ze 22)2
DRY (Zi,Zz)
RN j=1
10

---

Published in Transactions on Machine Learning Research (09/2023)
Obviously,
|A (Zi,Z>) L
DAN (Z,,2>) R. (Z1. Zz)
sup ZiEef 22E03
sup
Zi Er 22E05
AD (21, (22) D BRY (Z1, Zn)
+
sup
Z, EiE Z2El2
+ sup
DNN (Z,, L
ZIEC ZgEc2
7i T2 T3.
Lemma J.1, and equations (35) and (32) imply that there exists C > 0 such that
P[RN I RN S CNV/3= - 1.
Consequently, with probability going to 1, it happens that
Daw (Zi, Z2) < DAN (Zi, Z2) S DeN'S (2
Rx and trivially we have log N = o (CN'/) (CN)))) b) in Theorem 2.6 X,
Since, by assumption (15), log N o gives that T3 converges in probability to zero as N + co. Since L < co, this fact implics that
sup EIE ZEE0N
DRN Z2) 0p(1).
(42)
With respect to T2, we have
1
RN
(Zi
Z2, ,² 1A3 An1
T3
Z ElECT ZEEEN
sup
RN
f=1
Aj
Agj
AN
Aj Az1
DRy (Z, Z2) = Op(N O11101/10),
sup
Ø=1
Z ECN ZgEc
where last equality follows from (42), (32), (34), (35) and (17).
Finally, given Z1 E CN, Zz E C.N the Cauchy-Schward inequality and the fact that ||å = =Iløy|| = 1 imply
RN
(Z) Zzz 65 (Zi Zz Φ3)
|R. (2,, - D R. (Z Zz)
RN j=l
RN
(Z)
Zzz 63 أ,)ه
(Zi -Zz, 2, + +z))
RN
j=1
1
Ax
|6 - S 11 + +,11
<
Zzll
RN
1
Ry
11 di
< 2 Zzll
RN
2 7 Zzll H
Moreover, the application of (33), (34), (36) and (17) gives that HN = O00N--N11//0 which in turn is
equivalent to saying that there exists C' 0 such that Pilt, < CN CN-1/111 + 1. This and the reasoning
leading to (28) imply that to prove Ti P 0, it enough to show that for every C
22p [Z - Zoll > CN1/20] +0 as +0,
(43)
11

---

Published in Transactions on Machine Learning Research (09/2023)
where Z1 and Zy came from distributions P1 and F2, respectively.
To show (43), notice that 2 Z2 follows a Gaussian distribution whose mean function is 2 - and its
covariance is D = D2و + Eوz Let 7/ with j E N denote the ordered eigenvalues of 1 Consider a basis composed by the eigenfunctions of 1 we denote by (μ² - -2²2) the components of μ - م in this basis and {uj]eeN a sequence of i.i.d. real standard normal variables. Now, we have the following
1 Z2l 2
1/2 Uj + (n -12211)
j=1
20
Σ
(z5(1] 1) + 73 + (n12 -11 1122)3 + 2(p22 - -122) ///
j=1
20
Σ
(v(uz - 1) + 2( - 1 22) (2) + traco(E!2) + -11122 - 222
=
j=l
Note that K == trace(E!2) + |||221 -1.218 < 00. Thus,
8
7 -zzl² > CN1/10]
(vs(u} -1) + -1)+2(p² - 2 1/22//1 us) > CN1/100 K
= P
j=1
-1)> > 1/2 (CN1170- K)
S P
u5(u 1)
2
(1:2 11 2,7) uj > 1 (CN1//0011 -K)
+P
=: P P2.
(44)
-K) > 1 and from Lemma G.1, we have that Obviously, 1/2 (CN2/10- - K) + 6c. Thus, eventually 1/ (CN1/10- - K)
P lim P
Σ 2u(u3 -1)
1 2
(CN1/1O-K) 22exp (-3)
1 (CN1/1O-K)
(45)
For P2. first note that the real r.v. ΣØ (1:2 - -11222 ),²;) // un is centered normal, with variance equal to d
d
Σ (4 2r - 1. 3 I 7/ Σ (1 2h - -12203 S &71lpe? 21 - 2 for every d E N. Therefore,
1=1
1=1
d
P2
< d=000
lim P
Σ (112 - 1 1, /3 U3
1 1 (CNL/1O- x))
4
(CN1/10- (CN1/1O-K) - -k)
< P
1))
1/2 1 μ || 471
1 27/1 (4]1 - -122a1)2 (CN2/10-11 K) ))
exp
(46)
π
where last inequality comes from (25) because, eventually 1 1< (CN1/10 K) /447)/ /1 ()). Finally, (44), (45), and (46) give (43) and consequently, Ti L as N - xc.
12

---

Published in Transactions on Machine Learning Research (09/2023)
Appendix II: Additional Material
K Extension to non-Gaussian Processes
Obviously, non-Gaussian processes can also be mutually singular. In fact, Theorem 4.3 in Rao and Varadara- jan (1963) contains a sufficient condition for this property to be satisfied. This allows us to consider the possibility to extend previous results to cover non-Gaussian distributions. It is obvious that the proofs that we have developed can cover non-Gaussian distributions as long as they satisfy the due properties. In this subsection, we state the properties distribution should satisly in order the proofs can be extended. Let IP1 and P2 be two probabilities on the Hilbert space H. Here, Z will denote L2(0, 1] valueed random element with distribution Pi + x2P2 for some 11, 22 > 0 with क + 22 = 1.
The basic assumption is the existence of a covariance of Z. We will also consider assumptions A. and A.2 (see Section of the paper) and b E L2] Given a p.d. x matrix Ad and a d-dimensional subspace Va C L2(0, 1], we need to consider the d-dimensional random vector U, = (Aa) -1 4 b)a and the covariance matrix Sa= = =2/2//115/ where Σ/ is the covariance matrix of Zu and (Z b)a is the projection on Va of (Z b) with d N.
Let us write U₂ - U (u',.... udyT in the basis of the eigenvectors of Sy and let ......, ag denote the cigenvalues of Sa. Therefore, u,/a for I i S d are real, standardisod random variables which we ncod assume be i.dd. Similar properties must hold for the decomposition of Z its eigenfunction basis (also see Dai et al (2017)). We finally need two exponential inequalities those stated in (25) of Lemma G.I.
L Discussion on GP Clustering for the 'Location Only' Case
We have some idcas to fix the problem with the 'location only' case. Recall the notation used in Subsection 2.1 the paper. As stated there, the problem this case is that
D uu. v) = 2 (2 a vidl
Σ d
d
U 0;)2 P2 Oas d + oc.
A
Our idea is to replace the terms in the sum with some others going to 0 slowly enough (or, if possible, not converging to zero at all). To use this idea, our proposal is as follows:
D D (u, v) = (1))) v)a
1
with r T.
(47)
d
X
Here, II is the set of integers. In this article, we have studied the case when r= 1, i..,, Diae not strict requirement and we look into some possible scenarios below:
However, this
I E -1, then assumption A.2 in the main paper trivially gives that (u-- < (u,-
eventually for large i. Consequently, we have D
2a (u, v) f 0 as d co.
2 When E 3, the transformation D may be useful because 1/X, will start to take high values (recall assumption A.2) and this may lead to separation between the observations of corre- sponding to different clusters.
Keeping the viewpoint stated above in mind, we consider the transformation D (using r = 4 in (47)). Numerical results for the difference location only setting stated in Scction 4 are reported below. We have excluded Example II from this comparison because, as stated earlier, the two GPs have no differences in their means.
13

---

Published in Transactions on Machine Learning Research (09/2023)
Table L.1: Adjusted Rand distances for different GPs with differences only locations (with standard error in brackets).
CD
Ex.
k-means
spectral
mclust
funclust.
CT.
DHP
k- k-means
spectral
mclust
0.0001
0.0001
1.0000
0.0002
0.0001
0.0001
0. 0.0012
0.0814
0.0016
(0.0001)
(0.0001)
(0.0000)
(0.0001)
(0.0000)
(0.0001)
(0.0002)
(0.0072)
(0.0003)
III
0.0646
0.0960
1.0000
0. 0.0795
0.0945
0.1480
0.1965
0.1948
0.1649
(0.0015)
(0.0018)
(0.0000)
(0.0009)
(0.0045)
(0.0047)
(0.0025)
(0.0027)
(0.0017)
IV
0.1606
0.2111
1.0000
0.0318
0.1015
0.0134
0.1257
0.3777
0.1784
(0.0007)
(0.0008)
(0.0000)
(0.0003)
(0.0000)
(0.0004)
(0.0019)
(0.0085)
(0.0021)
As expected, the performance of k-means is quile good in Examples T and III. Both DHP and CL also perform quite well securing a first place in some cases. The proposed statistic D shows significant improvement (recall from part (b) of Theorem 2.2 that Lyg. = 0 for the earlier transformation D " and this is reflected in the numerical figures of Table L.1.
Clearly, there is scope of further work with the proposed transformation D Di for r E (2,3,...J both thooretically as well as numerically.
M Review of the paper by Delaigle et al (2019)
As stated, Delaigle et al (2019) is related to perfect clustering and it is the only paper on perfect clustering we are aware of. In this section, we analyze the relation between our proposal and this paper.
The proposal by Delaigle et al (2019) is based on finding finite-dimensional subspace in which the data are projected, and clustering is done by applying a modification of the k-means algorithm on those projections. A theoretical result related to perfect clustering is stated in Theorem of this paper. In the homoscedastic case, Delaigle et al (2019) gives an explicit expression of the subspace in which the data should be projected (see Theorem 2 of this paper).
The technique proposed in this paper has some advantages over our proposal in the sense that they can handle the homoscedastic (differences only in location) case. However, suffers from several limitations, the main being its difficulty to deal with more than two populations. As stated in Section 5.2 of Delaigle et al (2019), '...il is not clear how to ertend.. [their method] ...to more than two populations." 1 The only exception being the case in which one needs to assume that the data follow binary hierarchical structure, where the groups can be sequentially split into two. As consequence of this, Delaigle et al (2019) do not propose procedure to estimate the true number of clusters.
On the technical side, the theory of Delaigle et al (2019) has some limitations. It requires to arbitrarily fix p E N. Then, the data are projected on a p-dimensional subspace in which the clustering is be done. New issues appear in the way in which the subspace should be chosen as well as the clusters need to be constructed. According to Theorem 1 of this paper, the generators of the subspace must be chosen in finile set with cardinality on be chosen between those in finite set of Voronoi tessellations of TR with cardinality b, - as - co. Additionally, the result needs technical conditions like the existence of some E (0, 1) such that for every C it happens that bb exp(-Cro + 00 as n + 0o.
the sample size co. Moreover, the partition of the data set must
14

---

Published in Transactions on Machine Learning Research (09/2023)
N Full Numerical Results
Full results for two scenarios are given below.
Table N.1: Adjusted Rand distances for different GPs with differences in locations and scales (with standard error in brackets).
Ex.
CL1
CL2 0.0554
DHP1
DHP2
0.0239
0.8386
0.0818
(0.0007)
(0.0005)
(0.0064)
(0.0025)
0.5767
0.9967
0.5470
0.5149
(0.0045)
(0.0018)
(0.0047)
(0.0049)
0.2891
0.9962
0.4137
0.5613
(0.0000)
(0.0000)
(0.0054)
(0.0060)
0.1833
0.6660
0.1379
0.5211
(0.0000)
(0.0000)
(0.0033)
(0.0061)
Table N.2: Adjusted Rand distances for different GPs with differences only in scales (with standard error in brackets).
Ex.
CL1
CL2
DHP1
DHP2
0.8269
1,0017
0.9989
0.9966
(0.0000)
(0.0000)
(0.0005)
(0.0008)
0.9065
1.0019
1.0001
0.9999
(0.0007)
(0.0005)
(0.0002)
(0.0003)
0.9994
0.9998
0.9984
0.9967
(0.0000)
(0.0000)
(0.0004)
(0.0007)
IV
0.8464
0.9928
0.9994
0.9980
(0.0000)
(0.0000)
(0.0003)
(0.0006)
Full result for the location only case (using the transformation D is given next.
stated Section L Appendix II)
Table N.3: Adjusted Rand distances for different GPs with differences in locations (with standard error in brackets).
Ex.
C1
C2
DHP1
DHI2
0.0001
0.0001
0.9896
0.0001
(0.0000)
(0.0000)
(0.0001)
(0.0001)
III
0.0945
0.9975
0.1480
0.1623
(0.0010)
(0.0043)
(0.0017)
(0.0075)
IV
0,1015
0.9001
0.0134
0.1473
(0.0000)
(0.0000)
(0.0004)
(0.0039)
R codes for clustering methods available from this link: GP clustering.
15

---

Published in Transactions on Machine Learning Research (09/2023)
0 Pseudo-codes of Our Methods
Algorithm Clustering Algorithm
Data: Z...... ZN are observed at p time points (i.e.; p-dimensional data); Compute the sample variance-covarraacce matrix SN; Define dy min(p, N
D O(N2²)
Spectral decomposition of the dN dx matrix SN yields e Edy (eigen functions) and λ , Ady (eigen values); Fix < d < dx;
4 0(p")
for = : N do
for j= d do Compute 2ig = e,
end for end for for i= :1:NN: do for 1: do
D O(N3p)
d (z z
Compute Da(Zi, 71)
d k=1
A
end for end for
D O(N2d)
for i= 1 N do
for = N do
temp- = 0:
for l + i, do
Compute temp temp + [Da(Zz Zi, Z,) Da(Z,, 2))
end for Ta(Zz: Z2) temp;
end for end for Spectral decomposition of the x matrix I' yields ,I,...,Bu (eigen values);
b O(N2d) O(N3)
Use optishrink to get Ka;
Implement any suitable clustering algorithu on Px with Ka as input for the number of clusters; Obtain class labels for the observations, or a given data set.
Algorithm 2 Cross Validation Algorithm to Choose the Value Consider the data S [Z,...., Zn
for d = 2::0 dx do;
Fix B E N;
for 1- =:::A do
Split the set S into three disjoint subsets Sts, S25 and Sws Use Algorithm 1 with data points corresponding SIs to get cluster assignment for data points
corresponding to Sas (say, Ci);
Use Algorithm 1 with data points corresponding to S2, to get cluster assignment for data points corresponding to Sas (say, C2);
Compute distance between C1 and C2 (say, Ds) end for Average over BI...., Dg to get DCV; end for
D O(Bdx)
D O(dx)
Return dey = arg minggdedw DCV.
16


Reviews:
[{"review_id": "9kdH3oY1YP", "reviewer": "Reviewer_TJHS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The results are persuasive, and at a high level the paper is relatively straightforward to follow even without too much deep technical knowledge.\n- It is nice that what is a paper which (after some revision) could probably stand quite well only on its theoretical elements also contains some illustration for how it can be used practically, along with some reasonably satisfying performance results as well as code to reproduce the results\n\nWeaknesses:\n- My only major concern is that in trying to follow in detail, the notation could be made a lot clearer. This to the extent that it is not always clear that things are being correctly stated. I will detail some instances in the following section.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents some theoretical contributions relating to \"perfect clustering\" in the Gaussian Process (GP) mixture model framework for functional data. \n\nThe first set of results relate to the identifiability of the cluster \"labels\" of a fixed number of functions when viewed at finer and finer granularity via a sequence of associated finite dimensional distance metrics dependent on the covariance matrices of the projected mixture components, and the projection of the mixture itself. Because these distance functions depend only on the covariances their results do not directly allow for cluster identification in the scenario where only the mean functions of the components differ.\n The main result (from a clustering point of view) in this set states that the sequence of distances of the projections of two GPs tend in probability to distinct values depending on whether they are from two components with different covariance functions or not. Based on this they propose constructing a dissimilarity matrix in which the similarity between two processes is given as the average of the squared difference in distances between each of the two and each other \"realisation\" in the sample. Because of the previous result the dissimilarity will be zero (in probability) if they are from components with the same covariance and some non-zero value otherwise.\n\nThe second set of results relate to how this translates when the sequence of distance functions is based on sample covariance matrices in the increasing sample size context. The authors show that for (projection) dimension depending on the eigen-gaps of the sample covariance matrices they can essentially extend their distributional results extend to the sample context.\n\nSome discussion on how to realise their results (not all of which are directly constructive) can be utilised in practice, along with experimental results on simulated and benchmark data show some promising practical implications as well.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The statement of Theorem 2.1 is strange to me. Since (4) does not in any way depend on the function b, it should come before. Something like: Let {A_d} be ... for which \\lim_{d \\to \\infty} \\frac {||\\alpha_d||_\\infty} {d} = 0. Then if b is such that (2) and (3) we get our convergence result.\"\n- There is (as far as I can tell) a lack of consistency with when the covariance function (and its projections) is given its super-script Z (or Z_h depending on context) and when it is not. Am I missing some distinction here?\n- There is also inconsistency in when the covariance function in the superscript of D^{\\Sigma^Z_d}_d is given the subscript indicating it is actually a covariance function (otherwise it hasn't been defined when this supersctipt is not projected, even if there is no ambiguity in these cases).\n- In (8) why do you have both notation L^{hk}_S and also L^{hk}?\n- Statement (c) in Theorem 2.2 confuses me. In statements (a) and (b) you have had to assume the limits exist, whereas now suddenly these quantities suddenly exist and are also finite. Also, since h \\not = k why is there not an L^k_S which would be distinct from L^h_S and L^{hk}_S? There is no \"priority\" of h over k?\n- The definition of \\hat u^Z_j above (17) seems never to be used. Why is it included?\n- In the experiments, why is only the best result on the benchmark data sets given? All results for the simulations are given in the appendix and it would be interesting to know if the observation that the GMM applied on the \\Gamma matrices continued to be the most consistent over kmeans and spectral clustering.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "0Cg10OkG22", "reviewer": "Reviewer_ZAn5", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strength\n1. This paper is relevant as the problem of clustering and its theoretical analysis is still an important research topic.\n2. This paper carefully analyzed the problem of clustering based on the Hajek and Feldman property. As the authors argue in Introduction, to date, this approach has not been studied at the sufficient depth, and the contribution of this paper is novel.\n3. The proposed method shows good performance in experiments, which indicates that the contribution of this paper is not only theoretical but also potentially practical.\n\n## Weaknesses\n1. Overall this paper is hard to follow, and the presentation can be improved. One of the main reasons is that motivation and intuitive meaning has not been introduced when mathematical claims (Theorems) are stated. Thus it is difficult to understand why such theorems are required. In particular, in the context of clustering, how they are used in the task is often unclear.\n2. As I state in the \"requested changes\" section, there are several concerns regarding the validity of the assumption in the theoretical analysis. Therefore, in its current form, the quality and the significance of theoretical analysis is not convincing.\n3. Since experiments performed in this paper are not thorough, the quality of empirical evaluation is not high.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper theoretically analyzes data transformation under the assumption that data are generated from a mixture of Gaussian distributions. In particular, the task  of clustering is considered as a representative application and a practical clustering algorithm is proposed. Under some appropriate assumptions, the authors show that perfect clustering is realized when the sample size goes to infinity. Empirical evaluation on synthetic and real-world datasets shows that the proposal is superior to existing clustering methods in term of the adjusted Rand index.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Definition of $\\mathcal{C}\\_h$ in Eq.(6) is not clear to me. Now we are considering a single mixture distribution $\\mathbb{P}\\_{\\mathbf{Z}}$, and a sample should come from the distribution. Thus, although $\\mathbb{P}\\_1$ and $\\mathbb{P}\\_2$ are components of $\\mathbb{P}\\_{\\mathbf{Z}}$, the situation where a sample comes from either $\\mathbb{P}\\_1$ or $\\mathbb{P}\\_2$, which is considered in Eq.(6), does not make sense because this violates with the assumption of the mixture distribution with $0 < \\pi\\_1, \\pi\\_2 < 1$.\n2. Theorems in Section 2.1 seem to be based on the assumption of Eq.(6). Hence I think what the authors are treating is clustering of a sample from not a mixture distribution but two independent distributions. So I am not sure why Theorem 2.2 can be used in the example in Section 2.1.1.\n3. Around Eq.(6), $\\mathbf{Z}\\_1 ,\\dots, \\mathbf{Z}\\_N$ are used as elements of a sample, while $\\mathbf{Z}\\_1$ and $\\mathbf{Z}\\_2$ appear and are used as random processes (variables) in the sentence immediately before Theorem 2.2 as the authors use the notation $\\mathbb{P}\\_{\\mathbf{Z}\\_1}$ and $\\mathbb{P}\\_{\\mathbf{Z}\\_2}$. Are these $\\mathbf{Z}\\_1$ and $\\mathbf{Z}\\_2$ two elements of the sample, or some difference processes? If so, what are they in this context of clustering?\n4. Although the overview of the clustering process is given in Section 3.2, describing pseudo-code of the algorithm is also helpful for understanding.\n5. Please describe the computational complexity of the proposed algorithm.\n6. In Table 5.1, I guess $M$ means the sample size. Then please also write the number of features for each dataset. \n7. In experiments, it seems that what \"CD\" in text and tables means is not explicitly explained. So please explain it.\n8. All the real-world datasets in experiments are rather small scale. It would be interesting if larger scale, or higher-dimensional, datasets are also examined.", "broader_impact_concerns": "I do not have any concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9mbhpMX4hY", "reviewer": "Reviewer_BSFc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- This paper addressed an interesting problem: how to perfectly recover the component assignment for data generated from a mixture of Gaussian processes especially for the 2-component case through the lens of mutually singularity.\n\n- The authors provided solid theoretical foundations for both the proposed transformation and its application to Gaussian processes clustering, which lead to both estimation of the unknown number of clusters and the perfect separation of populations.\n\n- The experiments showcased the effectiveness of the proposed method comparing to different baseline methods.\n\nWeaknesses:\n- There is no proper definition of perfect clustering in the paper. As shown in [Priebe et. al. 2019], there is no universal definition of what a perfect clustering is. It seems that it refers to perfectly recover the component assignment if the data is generated from a mixture of Gaussian processes according to Definition 2.4 and Theorem 2.5. However, this definition cannot be directly applied to the J-mixture case with J > 2 since the proposed method produces more than J clusters as described in the first paragraph on page 9 which will not result an adjusted rand index approaching 1. Thus, it is unclear how to apply the proposed method to more than 2 clusters to achieve perfect clustering and the description presented is quite confusing. Also, there is no empirical validation on working with more than 2 clusters.\n\n- It would be better to specify the full set of assumptions/conditions for the main theorems. For example, the conditions used in Theorem 2.2 have been referenced multiple times in the paper, but they are not fully specified in the theorem. Additional assumptions were discussed in Section 2.1 as well. On the other hand, it is unclear if the choice of the $V_d$ as \"the d-dimensional subspaces generated by the d eigenfunctions associated with the d largest eigenvalues of $\\Sigma^{Z}$\" is critical to Theorem 2.2.\n\n- There lacks detailed full description of the proposed clustering algorithm. Even though the description in Section 3.2 appears quite simple, the actual implementation is much more involved as described in Section 4.2. It would be much better to have a single detailed description of the full implementation including ways to estimate the dimension d and the number of clusters J and highlight the differences between the theoretical derivations and empirical implementations.\n\n- Inconsistencies between theoretical derivations and empirical results. The derivations in Section 2 and 3 suggest that any clustering method can be applied on top of the transformed data to get perfect clustering and the authors also called out the k-means clustering there. However, there are large performance variations among different clustering methods applied to the transformed data especially in Example II and III as shown in Table N where k-means is usually the worst-performing one. This is particularly noteworthy since the theoretical analysis does not suggest such disparities. Unfortunately, the authors did not provide any discussion about why and how this happens in practice.\n\n- Unfair baseline comparisons and insufficient validations. First, the quantitative results reported in Section 4 and 5 are the best among 3 different clustering methods, which essentially does model selection on the test data. A more apples-to-apples comparison would be to compare k-means with CD-k-means and GMM with CD-mcluster, since the proposed method is focused on the transformed data not the clustering method on top. But when comparing k-means results in Table 4.1 with CD-k-means in Table N.1, we can see that vanilla k-means is actually performing better or comparable in three settings. Furthermore, the authors only considered the case with fixed number of samples and equal mixing weights of each component. It would be nicer to include results with different numbers of samples and varied mixing weights to showcase the asymptotic properties especially the relationship between N and d and handling of unequal mixing weights. Lastly, the proposed method seems to perform better in the \"difference in only scales\" case than the \"difference in location and scales\" case, even though the latter is an easier one. The authors might provide some discussions regarding this situation.\n\n\nPriebe et. al. 2019. On a two-truths phenomenon in spectral graph clustering. PNAS.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper first studied some transformation for Gaussian processes and its limit theorem through the lens of mutually singularity and then showed that it is possible to achieve perfect separations of mixtures of Gaussian processes through running common clustering algorithms like k-means on an appropriate modification of this transformation especially for the two components cases. This method performs better when comparing with different baseline methods on simulated and benchmark datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "On the presentation side,\n- Please consider rewriting the Abstract to not over-claim anything: the limitations of handling \"difference only in location\" case and if the proposed method only leads to perfect separation in 2-mixture. If the proposed method does lead to perfect separation in J-mixture with J>2, I would suggest having a dedicated section in the main text detailing the steps achieving that. The method used for 2-mixture does not seem to work with more than 2 components.\n- Please consider adding the problem formulation studied in this work along with the definition of perfect clustering in the Introduction, and including a Contributions subsection to make it easier for the readers to follow.\n- Please consider improving the presentation of technical part especially adding the full specification of the assumptions/conditions and providing more intuitive explanations for the Remarks. For example, it is not quite straightforward to see why b = 0 used in Remark 2.1.3 would satisfy all the conditions of Theorem 2.1.\n- Please consider improving the discussion in Remark 2.5.2. The argument of \"the structure ... will lead us to perfect clustering for every value of J $\\ge$ 2\" is not supported with any concrete method since the procedure described in the first paragraph on page 9 does not produce a perfect clustering. It also argues that \"the procedure described in Proposition 2.3 also works fine, with the limit equal to the rank of $\\Gamma$\", but the following paragraph mentioned that the rank does not generally coincides with the number of clusters.\n- Please consider adding the full description of the empirical implementations in Section 3.2 with all the details of estimating the dimension and the number of clusters. It would also be better to highlight the differences between the theoretical derivations and the empirical implementations.\n- Please consider correcting the discussion of Delaigle et al (2019) since the work extends their theoretical results to the case of more than two populations.\n\nOn the evaluation side,\n- Please consider including visualizations of typical curves from each group for each case to better illustrate the experimental settings and challenges.\n- Please consider addressing the inconsistencies between theoretical derivations and empirical results discussed above.\n- Please consider addressing the unfair baseline comparisons and insufficient validations discussed above.\n- Please consider adding empirical validation on working with more than 2 clusters depending on if the authors want to claim the perfect clustering in this scenario.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All three reviewers are leaning accept on the paper.  The revision in particular helped to clarify some of the concerns.  One reviewer still feels that the presentation of the paper could improve, and I would recommend doing a further pass to improve readability of the paper.   In particular, that reviewer stated "the paper would benefit from being proofread by a native English speaker as the grammar is at times awkward to read, not to mention sometimes incorrect."  I think in general the paper just needs a polishing pass---nothing specific has been noted by the reviewers but try to smooth it out. Otherwise, this looks like it's in good shape.
Decision: Accept (poster)