# Smooth Matérn Gaussian Random Fields: Euler Characteristic, Expected Number
and Height Distribution of Critical Points††thanks: Research partially
supported by NSF grant DMS-2220523 and the Simons Foundation Collaboration
Grant 854127.

Dan Cheng
Arizona State University

###### Abstract

This paper studies Gaussian random fields with Matérn covariance function with
smooth parameter $\nu>2$. Two cases of parameter spaces, the Euclidean space
and $N$-dimensional sphere, are considered. For such smooth Gaussian fields,
we have derived the explicit formulae for the expected Euler characteristic of
the excursion set, the expected number and height distribution of critical
points. The results are valuable for approximating the excursion probability
and computing p-values for peak inference.

Keywords: Gaussian random fields; Matérn; Smooth; Isotropic; Euler
characteristic; Height distribution; Critical points.

Mathematics Subject Classification: 60G15, 60G60, 62G32, 15B52.

##  1 Introduction

In recent years the Matérn class of covariance functions has gained widespread
popularity in spatial statistics [12, 14], mainly due to its great flexibility
for modelling the spatial covariance and hence dependent spatial data and
processes. Specifically, the Matérn covariance function is defined as

$\mathcal{M}(d)=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}\left(\frac{\sqrt{2\nu}d}{\ell}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}d}{\ell}\right),\quad
d\geq 0,$ (1.1)

where $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind.
This covariance function has three positive parameters $\sigma^{2}$, $\ell$
and $\nu$, with $\sigma^{2}$ controlling the variance, $1/\ell$ controlling
the spatial range, and $\nu$ controlling the smoothness. Let $X=\\{X(t),t\in
T\\}$ be a centered Gaussian random field living on an $N$-dimensional
parameter set $T$. We call $X$ a Matérn Gaussian random field if the
covariance satisfies

${\mathbb{E}}[X(t)X(s)]=\mathcal{M}(\|t-s\|),\quad t,s\in T.$

That is, $X$ is isotropic with covariance given by (1.1). There is rich
literature on Matérn Gaussian fields; however, only a few are focused on the
smooth case. This paper aims to bridging smooth Matérn Gaussian fields and
recent developments on the expected Euler characteristic of the excursion set
[1, 8], the expected number and height distribution of critical points [2, 5,
7, 6, 4].

Let $\\{X(t),t\in T\\}$ be a centered smooth Gaussian random field. Let
$A_{u}(X,T)=\\{t\in T:X(t)\geq u\\}$ be the excursion set of $X$ exceeding
level $u$ over the parameter set $T$ and denote by $\chi(A_{u}(X,T))$ its
Euler characteristic. It is shown in [1] that, the expected Euler
characteristic ${\mathbb{E}}[\chi(A_{u}(X,T))]$ is computable and can used to
approximate the excursion probability ${\mathbb{P}}[\sup_{t\in T}X(t)\geq u]$
for large $u$ such that the error is super-exponentially small. This is useful
in controlling the family-wise error in statistics [15].

The number of critical points of index $i$ of $X$ above $u$ is defined as

$\mu_{i}(X,u)=\\#\left\\{t\in D:X(t)\geq u,\nabla
X(t)=0,\text{index}(\nabla^{2}X(t))=i\right\\},\quad i=0,\ldots,N,$ (1.2)

where $D$ is an $N$-dimensional unit-area disc on $T$, $\nabla X(t)$ and
$\nabla^{2}X(t)$ are respectively the gradient and Hessian of $X$, and
$\text{index}(\nabla^{2}X(t))$ denotes the number of negative eigenvalues of
$\nabla^{2}X(t)$. Notice that we will focus on the expectation of
$\mu_{i}(X,u)$, which depends on the volume of $D$ (1 for a unit disc) due to
isotropy, thus we omit $D$ in the notation for simplicity. The expected number
of critical points ${\mathbb{E}}[\mu_{i}(X,u)]$ of smooth random fields is
important in statistics [6] and physics [2]. By default, let
$\mu_{i}(X)=\mu_{i}(X,-\infty)$ be the number of critical points of index $i$
of $X$ over $D$.

The height distribution of a critical value of index $i$ of $X$ at $t_{0}$ is
defined as

$F_{i}(u):=\lim_{\varepsilon\to 0}{\mathbb{P}}\left\\{X(t_{0})>u\ |\
\exists\text{ a critical point of index $i$ of }X(t)\text{ in
}B(t_{0},\varepsilon)\right\\},$ (1.3)

where $B(t_{0},\varepsilon)$ is the ball of radius $\varepsilon$ centered at
$t_{0}$. It has been found to be an important tool for computing p-values in
peak detection and thresholding problems in statistics [7, 3, 4] and
neuroimaging [9]. Notice that, due to isotropy, $F_{i}(u)$ does not depend on
the location $t_{0}$. It is shown in [5] that,

$F_{i}(u)=\frac{{\mathbb{E}}[\mu_{i}(X,u)]}{{\mathbb{E}}[\mu_{i}(X)]}.$ (1.4)

Therefore, $F_{i}(u)$ can be obtained immediately once the form of
${\mathbb{E}}[\mu_{i}(X,u)]$ is known.

##  2 Smooth Matérn Gaussian Random Fields on Euclidean Space

It is well known (see for example [14]) that the Gaussian fields with
covariance function given by (1.1) are $m$ times differentiable if and only if
$\nu>m$. To study the expected Euler characteristic and critical points, it
requires the Gaussian fields to be twice differentiable [1, 5]; thus,
throughout this paper, we assume $\nu>2$.

To characterize the variances of $X^{\prime}$ and $X^{\prime\prime}$,
following the notation in [5], we introduce a real function
$\rho(d)=\mathcal{M}(d^{2})$ and have the following results.

######  Proposition 2.1.

Let $\rho(d)=\mathcal{M}(d^{2})$. Suppose $\nu>2$. Then

$\begin{split}\rho^{\prime}(0)=-\frac{\sigma^{2}\nu}{2(\nu-1)\ell^{2}}\quad\text{\rm
and
}\quad\rho^{\prime\prime}(0)=\frac{\sigma^{2}\nu^{2}}{4(\nu-1)(\nu-2)\ell^{4}}.\end{split}$
(2.1)

In particular,

$\begin{split}\kappa:=-\frac{\rho^{\prime}(0)}{\sqrt{\rho^{\prime\prime}(0)}}\bigg{\rvert}_{\sigma=1}=\sqrt{\frac{\nu-2}{\nu-1}}\quad\text{\rm
and
}\quad\eta:=\frac{\sqrt{-\rho^{\prime}(0)}}{\sqrt{\rho^{\prime\prime}(0)}}=\sqrt{\frac{2(\nu-2)}{\nu}}\ell.\end{split}$
(2.2)

###### Proof.

It follows from the property of Bessel functions (see page 502 in [13]) that,
for non-integer $\nu$,

$\begin{split}K_{\nu}(r)=\frac{2^{\nu-1}\Gamma(\nu)}{r^{\nu}}\sum_{j=0}^{\infty}\frac{(r^{2}/4)^{j}}{j!(1-\nu)_{j}}+\frac{r^{\nu}\Gamma(-\nu)}{2^{\nu+1}}\sum_{j=0}^{\infty}\frac{(r^{2}/4)^{j}}{j!(1+\nu)_{j}},\end{split}$
(2.3)

where $(x)_{j}=x(x+1)\cdots(x+j-1)$ with $(x)_{0}=1$ is the Pochhammer symbol.
For $\nu>2$, organizing the right side of (2.3), we obtain the following
expansion as $r\to 0$,

$\begin{split}\frac{r^{\nu}}{2^{\nu-1}\Gamma(\nu)}K_{\nu}(r)=1-\frac{r^{2}}{4(\nu-1)}+\frac{r^{4}}{32(\nu-1)(\nu-2)}+o(r^{4}).\end{split}$
(2.4)

On the other hand, by page 502 in [13], if $\nu=n$ is an integer, then

$\begin{split}K_{n}(r)=&\frac{2^{n-1}}{r^{n}}\sum_{j=0}^{n-1}\frac{(-1)^{j}(n-j-1)!(r^{2}/4)^{j}}{j!}\\\
&+\frac{(-1)^{n}r^{n}}{2^{n}}\sum_{j=0}^{\infty}\left[\frac{\psi(1+j)}{2}+\frac{\psi(1+n+j)}{2}-\log\left(\frac{r}{2}\right)\right]\frac{(r^{2}/4)^{j}}{j!(n+j)!},\end{split}$

where $\psi(\cdot)$ is the digamma function defined by
$\psi(z)=\frac{d}{dz}\log\Gamma(z)=\Gamma^{\prime}(z)/\Gamma(z)$. This implies
that, for $n\geq 3$, we have the following expansion as $r\to 0$,

$\begin{split}\frac{r^{n}}{2^{n-1}\Gamma(n)}K_{n}(r)=1-\frac{r^{2}}{4(n-1)}+\frac{r^{4}}{32(n-1)(n-2)}+o(r^{4}).\end{split}$
(2.5)

Combining (2.4) and (2.5), we obtain that, for $\nu>2$,

$\begin{split}\mathcal{M}(d)&=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}\left(\frac{\sqrt{2\nu}d}{\ell}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}d}{\ell}\right)\\\
&=\sigma^{2}\left[1-\frac{\nu}{2(\nu-1)}\left(\frac{d}{\ell}\right)^{2}+\frac{\nu^{2}}{8(\nu-1)(\nu-2)}\left(\frac{d}{\ell}\right)^{4}\right]+o(d^{4}),\quad
d\to 0.\end{split}$ (2.6)

In other words, we have

$\begin{split}\rho(d)=\sigma^{2}\left[1-\frac{\nu}{2(\nu-1)\ell^{2}}d+\frac{\nu^{2}}{8(\nu-1)(\nu-2)\ell^{4}}d^{2}\right]+o(d^{2}),\quad
d\to 0.\end{split}$

This second-order Taylor expansion for $\rho(d)$ around $d=0$ implies (2.1)
and hence (2.2). ∎

###  2.1 The expected Euler characteristic of excursion set

Denote by $\phi(x)=(2\pi)^{-1/2}e^{-x^{2}/2}$ and
$\Psi(x)=\int_{x}^{\infty}\phi(x)dx$ the density and tail probability of the
standard normal distribution, respectively. Let $H_{j}(x)$ be the Hermite
polynomial of order $j$, i.e.,

$H_{j}(x)=(-1)^{j}e^{x^{2}/2}\frac{d^{j}}{dx^{j}}\big{(}e^{-x^{2}/2}\big{)},\quad
j\geq 0.$

######  Theorem 2.2.

Let $\\{X(t),t\in T\\}$ be a centered Gaussian random field with Matérn
covariance function (1.1), where $T\subset{\mathbb{R}}^{N}$ is an
$N$-dimensional piecewise smooth set. Suppose $\nu>2$. Then the expected Euler
characteristic is given by

$\begin{split}{\mathbb{E}}[\chi(A_{u}(X,T))]=\sum_{j=0}^{N}\frac{\nu^{j/2}}{(\nu-1)^{j/2}\ell^{j}}\mathcal{L}_{j}(T)\xi_{j}\left(\frac{u}{\sigma}\right),\end{split}$
(2.7)

where $\mathcal{L}_{j}(T)$ are the Lipschitz-Killing curvatures (cf. (10.7.3)
in [1]) of $T$ and

$\begin{split}\xi_{0}(x)=\Psi(x),\quad\xi_{j}(x)=(2\pi)^{-j/2}H_{j-1}(x)\phi(x),\quad
j\geq 1.\end{split}$ (2.8)

###### Proof.

It follows from formula (5.5.5) in [1] or Lemma 3.2 in [6] that ${\rm
Var}[\nabla X(t)]=-2\rho^{\prime}(0)I_{N}$, where $I_{N}$ is the $N\times N$
identity matrix. By (2.1), we obtain

${\rm Var}[\nabla X(t)/\sigma]=\left(\frac{\nu}{(\nu-1)\ell^{2}}\right)I_{N}.$
(2.9)

Applying Theorem 12.4.1 in [1] to the standardized Gaussian field $X/\sigma$
yields

$\begin{split}{\mathbb{E}}[\chi(A_{u}(X,T))]={\mathbb{E}}[\chi(A_{u/\sigma}(X/\sigma,T))]=\sum_{j=0}^{N}\mathcal{L}_{j}^{X/\sigma}(T)\xi_{j}\left(\frac{u}{\sigma}\right),\end{split}$
(2.10)

where $\mathcal{L}_{j}^{X/\sigma}(T)$ are the Lipschitz-Killing curvatures of
$T$ calculated with respect to the metric (cf. (12.2.2) in [1]) induced by
$X/\sigma$. Due to (2.9), following the arguments on page 423 in [1], we see
that the metric induced by $X/\sigma$ makes a new inner product for
$t,s\in{\mathbb{R}}^{N}$ given by
$\frac{\nu}{(\nu-1)\ell^{2}}{\langle}t,s\rangle$, where
${\langle}\cdot,\cdot\rangle$ is the simple Euclidean inner product, implying
that

$\mathcal{L}_{j}^{X/\sigma}(T)=\left(\frac{\nu}{(\nu-1)\ell^{2}}\right)^{j/2}\mathcal{L}_{j}(T),\quad
j\geq 0.$

plugging this into (2.10) yields the desired result (2.7). ∎

The formula (2.7) shows that the expected Euler characteristic is computable
with a relatively simple form, mainly due to the isotropy of Matérn Gaussian
fields. The Lipschitz-Killing curvatures $\mathcal{L}_{j}(T)$ depend on the
geometry of $T$. We show below an example for the case when $T$ is a cube.

Example 2.3  Let $T=[0,b]^{N}$ be an $N$-dimensional cube in
${\mathbb{R}}^{N}$. Then, by (10.7.4) in [1], $\mathcal{L}_{j}(T)={N\choose
j}b^{j}$, which implies

$\begin{split}{\mathbb{E}}[\chi(A_{u}(X,T))]=\phi\left(\frac{u}{\sigma}\right)\sum_{j=1}^{N}\frac{{N\choose
j}b^{j}\nu^{j/2}}{(2\pi)^{j/2}(\nu-1)^{j/2}\ell^{j}}H_{j-1}\left(\frac{u}{\sigma}\right)+\Psi\left(\frac{u}{\sigma}\right).\end{split}$

###  2.2 Expected number and height distribution of critical points

It is introduced in [6] that, an $N\times N$ random matrix $M=(M_{ij})_{1\leq
i,j\leq N}$ is called _Gaussian Orthogonally Invariant_ (GOI) with _covariance
parameter_ $c$, denoted by ${\rm GOI}(c)$, if it is symmetric and all entries
are centered Gaussian variables such that

${\mathbb{E}}[M_{ij}M_{kl}]=\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+c\delta_{ij}\delta_{kl},$
(2.11)

where $\delta_{ij}$ is the Kronecker delta function. We see that ${\rm
GOI}(c)$ becomes a GOE matrix if $c=0$. In particular, the density of the
ordered eigenvalues $\lambda_{1}\leq\ldots\leq\lambda_{N}$ of ${\rm GOI}(c)$
is

$\begin{split}f_{c}(\lambda_{1},\ldots,\lambda_{N})&=\frac{1}{K_{N}\sqrt{1+Nc}}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{N}\lambda_{i}^{2}+\frac{c}{2(1+Nc)}\left(\sum_{i=1}^{N}\lambda_{i}\right)^{2}\right\\}\\\
&\quad\times\prod_{1\leq i<j\leq
N}|\lambda_{i}-\lambda_{j}|\mathbbm{1}_{\\{\lambda_{1}\leq\ldots\leq\lambda_{N}\\}},\end{split}$
(2.12)

where $K_{N}=2^{N/2}\prod_{i=1}^{N}\Gamma\left(\frac{i}{2}\right)$ and
$c>-1/N$. We use the notation ${\mathbb{E}}_{{\rm GOI}(c)}^{N}$ to represent
the expectation under the ${\rm GOI}(c)$ density (2.12), i.e., for a
measurable function $g$,

${\mathbb{E}}_{{\rm
GOI}(c)}^{N}[g(\lambda_{1},\ldots,\lambda_{N})]=\int_{{\mathbb{R}}^{N}}g(\lambda_{1},\ldots,\lambda_{N})f_{c}(\lambda_{1},\ldots,\lambda_{N})d\lambda_{1}\cdots
d\lambda_{N}.$ (2.13)

######  Theorem 2.4.

Let $\\{X(t),t\in T\\}$ be a centered Gaussian random field with Matérn
covariance function (1.1), where $T\subset{\mathbb{R}}^{N}$ is an
$N$-dimensional set. Suppose $\nu>2$. Then for $i=0,\ldots,N$,

$\begin{split}{\mathbb{E}}[\mu_{i}(X)]&=\frac{2^{N/2}}{\pi^{N/2}\eta^{N}}{\mathbb{E}}_{{\rm
GOI}(1/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{i}<0<\lambda_{i+1}\\}}\right],\\\
{\mathbb{E}}[\mu_{i}(X,u)]&=\frac{2^{N/2}}{\pi^{N/2}\eta^{N}}\int_{u/\sigma}^{\infty}\phi(x){\mathbb{E}}_{{\rm
GOI}((1-\kappa^{2})/2)}^{N}\left[\prod_{j=1}^{N}\big{|}\lambda_{j}-\kappa
x/\sqrt{2}\big{|}\mathbbm{1}_{\\{\lambda_{i}<\kappa
x/\sqrt{2}<\lambda_{i+1}\\}}\right]dx,\\\
F_{i}(u)&=\frac{\int_{u/\sigma}^{\infty}\phi(x){\mathbb{E}}_{{\rm
GOI}((1-\kappa^{2})/2)}^{N}\left[\prod_{j=1}^{N}\big{|}\lambda_{j}-\kappa
x/\sqrt{2}\big{|}\mathbbm{1}_{\\{\lambda_{i}<\kappa
x/\sqrt{2}<\lambda_{i+1}\\}}\right]dx}{{\mathbb{E}}_{{\rm
GOI}(1/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{i}<0<\lambda_{i+1}\\}}\right]},\end{split}$
(2.14)

where $\kappa$ and $\eta$ are given in (2.2), ${\mathbb{E}}_{{\rm
GOI}(c)}^{N}$ is defined in (2.12) and $\lambda_{0}$ and $\lambda_{N+1}$ are
regarded respectively as $-\infty$ and $\infty$ for consistency.

###### Proof.

Note that, by the definition (1.2), we have $\mu_{i}(X)=\mu_{i}(X/\sigma)$ and
$\mu_{i}(X,u)=\mu_{i}(X/\sigma,u/\sigma)$. Applying Theorem 3.5 in [6] to $X$
and the standardized field $X/\sigma$ respectively, we obtain the first and
second lines in (2.14). Finally, the last line in (2.14) follows directly from
(1.4). ∎

##  3 Smooth Matérn Gaussian Random Fields on Spheres

The applications in geoscience, astronomy and environmental sciences have
stimulated recent rapid development in statistics of random fields on spheres.
It has been shown in [10] that, many of the commonly used covariance functions
on Euclidean spaces are valid on spheres when Euclidean distance is replaced
by the spherical distance (great circle distance) on a sphere. However, the
Matérn class in (1.1) is positive definite with the spherical distance only if
$\nu\leq 1/2$.

Let ${\mathbb{S}}^{N}$ be an $N$-dimensional unit sphere. Let $\|\cdot\|$ and
${\langle}\cdot,\cdot\rangle$ be the Euclidean distance and inner product in
${\mathbb{R}}^{N+1}$, respectively. As shown in [16], one can apply the
identity

$\|x-y\|=2\sin\left(\frac{\theta(x,y)}{2}\right),\quad\forall
x,\,y\in\mathbb{S}^{N}\subset{\mathbb{R}}^{N+1},$

where $\theta(x,y)=\arccos{\langle}x,y\rangle\in[0,\pi]$ denotes the spherical
distance on $\mathbb{S}^{N}$, to construct covariance functions on spheres. In
particular, by (1.1), we define the Matérn covariance function on
${\mathbb{S}}^{N}$ as

$\mathcal{M}\left(2\sin(\theta(x,y)/2)\right)=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}\left(\frac{2\sqrt{2\nu}\sin(\theta(x,y)/2)}{\ell}\right)^{\nu}K_{\nu}\left(\frac{2\sqrt{2\nu}\sin(\theta(x,y)/2)}{\ell}\right).$
(3.1)

Let $\\{X(x),x\in{\mathbb{S}}^{N}\\}$ be a centered Gaussian random field with
Matérn covariance function (3.1). Then we have

${\mathbb{E}}[X(x)X(y)]=\mathcal{M}\left(2\sin(\theta(x,y)/2)\right),\quad\forall
x,\,y\in\mathbb{S}^{N}.$ (3.2)

By Corollary 1 in [11], one still has the property that the Gaussian field $X$
is twice differentiable if and only if $\nu>2$. Note that

$\sin(\theta(x,y)/2)=\sqrt{\frac{1-\cos\theta(x,y)}{2}}=\sqrt{\frac{1-{\langle}x,y\rangle}{2}}.$

Following the notation in [7], we write the covariance in (3.2) as

$\begin{split}C({\langle}x,y\rangle)&:={\mathbb{E}}[X(x)X(y)]=\mathcal{M}\left(\sqrt{2(1-{\langle}x,y\rangle)}\right)\\\
&=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}\left(\frac{\sqrt{2\nu}\sqrt{2(1-{\langle}x,y\rangle)}}{\ell}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}\sqrt{2(1-{\langle}x,y\rangle)}}{\ell}\right),\
\forall x,\,y\in\mathbb{S}^{N},\end{split}$ (3.3)

and derive the following results.

######  Proposition 3.1.

Let $\\{X(x),x\in{\mathbb{S}}^{N}\\}$ be a centered Gaussian random field with
Matérn covariance function (3.3). Suppose $\nu>2$. Let
$C(p):=C({\langle}x,y\rangle)$, $p\in[-1,1]$. Then

$\begin{split}C^{\prime}(1)=\frac{\sigma^{2}\nu}{(\nu-1)\ell^{2}}\quad\text{\rm
and }\quad
C^{\prime\prime}(1)=\frac{\sigma^{2}\nu^{2}}{(\nu-1)(\nu-2)\ell^{4}}.\end{split}$
(3.4)

In particular,

$\begin{split}\tilde{\kappa}:=\frac{C^{\prime}(1)}{\sqrt{C^{\prime\prime}(1)}}\bigg{\rvert}_{\sigma=1}=\sqrt{\frac{\nu-2}{\nu-1}}\quad\text{\rm
and
}\quad\tilde{\eta}:=\frac{\sqrt{C^{\prime}(1)}}{\sqrt{C^{\prime\prime}(1)}}=\frac{\sqrt{\nu-2}}{\sqrt{\nu}}\ell.\end{split}$
(3.5)

###### Proof.

By (2.6), we obtain that, as $p\to 1$,

$\begin{split}C(p)&=\mathcal{M}\left(\sqrt{2(1-p)}\right)\\\
&=\sigma^{2}\left[1-\frac{\nu}{(\nu-1)\ell^{2}}(1-p)+\frac{\nu^{2}}{2(\nu-1)(\nu-2)\ell^{4}}(1-p)^{2}\right]+o\left((1-p)^{2}\right)\\\
&=\sigma^{2}\left[1+\frac{\nu}{(\nu-1)\ell^{2}}(p-1)+\frac{\nu^{2}}{2(\nu-1)(\nu-2)\ell^{4}}(p-1)^{2}\right]+o\left((p-1)^{2}\right).\end{split}$
(3.6)

This second-order Taylor expansion for $C(p)$ around $p=1$ implies (3.4) and
hence (3.5). ∎

Let $\omega_{j}=\frac{2\pi^{(j+1)/2}}{\Gamma((j+1)/{2})}$ be the spherical
area of the $j$-dimensional unit sphere $\mathbb{S}^{j}$. We have the
following results on the expected Euler characteristic of the excursion set
$A_{u}(X,{\mathbb{S}}^{N})=\\{x\in{\mathbb{S}}^{N}:X(x)\geq u\\}$, as well as
the expected number and height distribution of critical points.

######  Theorem 3.2.

Let $\\{X(x),x\in{\mathbb{S}}^{N}\\}$ be a centered Gaussian random field with
Matérn covariance function (3.3). Suppose $\nu>2$. Then the expected Euler
characteristic is

$\begin{split}{\mathbb{E}}[\chi(A_{u}(X,{\mathbb{S}}^{N}))]=\sum_{j=0}^{N}\frac{\nu^{j/2}}{(\nu-1)^{j/2}\ell^{j}}\mathcal{L}_{j}(\mathbb{S}^{N})\xi_{j}\left(\frac{u}{\sigma}\right),\end{split}$
(3.7)

where $\xi_{j}(\cdot)$ are given in (2.8) and

$\begin{split}\mathcal{L}_{j}(\mathbb{S}^{N})=\left\\{\begin{array}[]{l
l}2\binom{N}{j}\frac{\omega_{N}}{\omega_{N-j}}&\quad\text{if $N-j$ is
even,}\\\ 0&\quad\text{otherwise},\end{array}\right.\quad
j=0,1,\ldots,N,\end{split}$ (3.8)

are the Lipschitz-Killing curvatures of $\mathbb{S}^{N}$.

###### Proof.

Applying Lemma 3.5 in [8] to the standardized Gaussian field $X/\sigma$ yields

$\begin{split}{\mathbb{E}}[\chi(A_{u}(X,{\mathbb{S}}^{N}))]&={\mathbb{E}}[\chi(A_{u/\sigma}(X/\sigma,{\mathbb{S}}^{N}))]=\sum_{j=0}^{N}\left(\frac{C^{\prime}(1)}{\sigma^{2}}\right)^{j/2}\mathcal{L}_{j}(\mathbb{S}^{N})\xi_{j}\left(\frac{u}{\sigma}\right)\\\
&=\sum_{j=0}^{N}\left(\frac{\nu}{(\nu-1)\ell^{2}}\right)^{j/2}\mathcal{L}_{j}(\mathbb{S}^{N})\xi_{j}\left(\frac{u}{\sigma}\right),\end{split}$

where $C^{\prime}(1)$ is given in (3.4). ∎

######  Theorem 3.3.

Let $\\{X(x),x\in{\mathbb{S}}^{N}\\}$ be a centered Gaussian random field with
Matérn covariance function (3.3). Suppose $\nu>2$. Then for $i=0,\ldots,N$,

$\begin{split}{\mathbb{E}}[\mu_{i}(X)]&=\frac{1}{\pi^{N/2}\tilde{\eta}^{N}}{\mathbb{E}}_{{\rm
GOI}((1+\tilde{\eta}^{2})/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{i}<0<\lambda_{i+1}\\}}\right],\\\
{\mathbb{E}}[\mu_{i}(X,u)]&=\frac{1}{\pi^{N/2}\tilde{\eta}^{N}}\int_{u/\sigma}^{\infty}\phi(x){\mathbb{E}}_{{\rm
GOI}((1+\tilde{\eta}^{2}-\tilde{\kappa}^{2})/2)}^{N}\left[\prod_{j=1}^{N}\big{|}\lambda_{j}-\tilde{\kappa}x/\sqrt{2}\big{|}\mathbbm{1}_{\\{\lambda_{i}<\tilde{\kappa}x/\sqrt{2}<\lambda_{i+1}\\}}\right]dx,\\\
F_{i}(u)&=\frac{\int_{u/\sigma}^{\infty}\phi(x){\mathbb{E}}_{{\rm
GOI}((1+\tilde{\eta}^{2}-\tilde{\kappa}^{2})/2)}^{N}\left[\prod_{j=1}^{N}\big{|}\lambda_{j}-\tilde{\kappa}x/\sqrt{2}\big{|}\mathbbm{1}_{\\{\lambda_{i}<\tilde{\kappa}x/\sqrt{2}<\lambda_{i+1}\\}}\right]dx}{{\mathbb{E}}_{{\rm
GOI}((1+\tilde{\eta}^{2})/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{i}<0<\lambda_{i+1}\\}}\right]},\end{split}$
(3.9)

where $\tilde{\kappa}$ and $\tilde{\eta}$ are given in (3.5),
${\mathbb{E}}_{{\rm GOI}(c)}^{N}$ is defined in (2.12) and $\lambda_{0}$ and
$\lambda_{N+1}$ are regarded respectively as $-\infty$ and $\infty$ for
consistency.

###### Proof.

Applying Theorem 4.4 in [6] to $X$ and the standardized field $X/\sigma$
respectively, we obtain the first and second lines in (3.9). The last line in
(3.9) is a direct consequence of (1.4). ∎

## References

  * [1] Adler, R. J. and Taylor, J. E. (2007). Random Fields and Geometry. Springer, New York.
  * [2] Bardeen, J. M., Bond, J. R., Kaiser, N. and Szalay A. S. (1985). The statistics of peaks of Gaussian random fields. Astrophys. J., 304, 15–61.
  * [3] Cheng, D., Cammarota, V., Fantaye, Y., Marinucci, D. and Schwartzman, A. (2020). Multiple testing of local maxima for detection of peaks on the (celestial) sphere. Bernoulli, 26, 31–60.
  * [4] Cheng, D., He, Z. and Schwartzman, A. (2020). Multiple testing of local extrema for detection of change points. Electron. J. Statist., 14, 3705–3729.
  * [5] Cheng, D. and Schwartzman, A. (2015). Distribution of the height of local maxima of Gaussian random fields. Extremes, 18, 213–240.
  * [6] Cheng, D. and Schwartzman, A. (2018). Expected number and height distribution of critical points of smooth isotropic Gaussian random fields. Bernoulli, 24, 3422–3446.
  * [7] Cheng, D. and Schwartzman, A. (2017). Multiple testing of local maxima for detection of peaks in random fields. Ann. Stat., 45, 529–556.
  * [8] Cheng, D. and Xiao, Y. (2016), Excursion probability of Gaussian random fields on sphere. Bernoulli, 22, 1113–1130.
  * [9] Chumbley, J. R and Friston, K. J. (2009). False discovery rate revisited: FDR and topological inference using Gaussian random fields. _Neuroimage_ , 44, 62–70.
  * [10] Gneiting, T. (2013), Strictly and non-strictly positive definite functions on spheres. Bernoulli 19, 1327–1349.
  * [11] Guinness, J. and Fuentes, M. (2016), Isotropic covariance functions on spheres: Some properties and modeling considerations. Journal of Multivariate Analysis 143, 143–152.
  * [12] Guttorp, P. and Gneiting, T. (2006). Studies in the history of probability and statistics XLIX: On the Matérn correlation family. Biometrika. 93, 989–995.
  * [13] Spanier, J. and Oldham, K. (1987). An Atlas of Functions. Springer.
  * [14] Stein, M. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer.
  * [15] Taylor, J. E. and Worsley, K. J. (2007). Detecting sparse signals in random fields, with an application to brain mapping. _J. Am. Statist. Assoc._ , 102, 913–928.
  * [16] Yadrenko, A. M. (1983). Spectral Theory of Random Fields. Optimization Software, New York.

> Dan Cheng
> School of Mathematical and Statistical Sciences
> Arizona State University
> 900 S. Palm Walk
> Tempe, AZ 85281, U.S.A.
> E-mail<EMAIL_ADDRESS>


# LSEC: Large-scale spectral ensemble clustering

Hongmin Li, Xiucai Ye, Akira Imakura and Tetsuya Sakurai

(Department of Computer Science
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>{yexiucai, imakura,
<EMAIL_ADDRESS>
)

###### Abstract

Ensemble clustering is a fundamental problem in the machine learning field,
combining multiple base clusterings into a better clustering result. However,
most of the existing methods are unsuitable for large-scale ensemble
clustering tasks due to the efficiency bottleneck. In this paper, we propose a
large-scale spectral ensemble clustering (LSEC) method to strike a good
balance between efficiency and effectiveness. In LSEC, a large-scale spectral
clustering based efficient ensemble generation framework is designed to
generate various base clusterings within a low computational complexity. Then
all based clustering are combined through a bipartite graph partition based
consensus function into a better consensus clustering result. The LSEC method
achieves a lower computational complexity than most existing ensemble
clustering methods. Experiments conducted on ten large-scale datasets show the
efficiency and effectiveness of the LSEC method. The MATLAB code of the
proposed method and experimental datasets are available at
https://github.com/Li-Hongmin/MyPaperWithCode.

##  1 Introduction

Ensemble clustering, also known as consensus clustering, is a classic problem
in machine learning field, aiming to combine multiple base clustering into a
better and more consensus clustering [24, 18, 6, 14, 28, 22, 32, 11, 12, 21,
20, 8, 10, 37, 38]. Due to its good performance, ensemble clustering has a
pivotal role in many research areas, such as community detection [25] and
bioinformatics [15, 29].

There are two critical steps in ensemble clustering: ensemble generation and
consensus function. Ensemble generation aims to generate multiple base
clusterings on the same datasets. In the early stage, $k$-means based ensemble
generation methods [14, 27, 20] are widely used. Recently, spectral clustering
based ensemble generation methods [13, 16] have received attention for its
high performance. On the other hand, the consensus function is used to
integrate multiple base clusterings into a consensus one. We can roughly
categorize ensemble clustering according to the consensus function into two
categories: the co-association matrix based methods and the graph partitioning
based methods.

The co-association matrix based ensemble clustering method [6, 14, 31, 30] is
one of the most widely used ensemble clustering strategies. A typical example
is the evidence accumulation clustering method [6], which counts the frequency
of the pair-wise co-occurrence of the same cluster between a pair data points
according to base clusterings. After treating the co-association matrix as a
similarity matrix, the hierarchical agglomerative clustering algorithm is
applied to obtain the consensus clustering. Iam-On et al. [14] extend the EAC
method by constructing the co-association matrix based on the common
neighborhood information between clusters. Tao et al. [26] propose a robust
spectral ensemble clustering method to learns a robust representation for the
co-association matrix by capturing the noises and conduct spectral clustering
to obtain consensus clustering. Huang et al. [12] also enhance the co-
association matrix based on similarity mapping from the cluster-level to the
object-level and achieve ensemble clustering via fast propagation of cluster-
wise similarities. However, the co-occurrence matrix based methods often lead
to high computational cost, which has become a bottleneck for large-scale
clustering tasks. Therefore, most co-association matrix based methods can work
well in small-scale datasets but hardly complete large-scale clustering tasks
in an acceptable time.

Graph partitioning based ensemble clustering methods [24, 5, 9, 16] aim to
transform the ensemble clustering problem into a graph partitioning problem to
find the consensus clustering. Strehl and Ghosh [24] construct a hypergraph
representation by exploring base clusterings and propose three graph
partitioning based ensemble clustering methods. Huang et al. [9] develop a
sparse graph with a small number of probably reliable links from base
clusterings and find the consensus clustering based on probability trajectory
analysis. Li et al. [16] apply spectral clustering method as base clusterings
and take the graph Laplacian matrices of base clusterings as input, then learn
a consensus representation by optimizing the graph Laplacians of consensus
clustering and base clusterings simultaneously, finally conduct spectral
clustering to obtain consensus clustering. Although graph partitioning based
methods have successfully improved clustering quality, they still have
limitations regarding large-scale datasets.

Recently, a few studies have made progress in the application of large-scale
data for ensemble clustering. Wu et al. [32] propose a $k$-means based
consensus clustering (KCC) method, which applies the $k$-means method on a
contingency matrix from base clusterings to obtain the consensus clustering
result efficiently. Liu et al. [20] transform the spectral clustering of the
co-association matrix into a weighted $k$-means method and prove two
approaches are equivalent, which achieve high efficiency for ensemble spectral
clustering. Huang et al. [13] point out the efficient bottleneck of $k$-means
based ensemble generation and apply a large-scale spectral clustering method
to fast product the base clusterings, then conduct bipartite graph
partitioning to obtain the consensus clustering. Although these studies have
achieved success in their respective fields, the large-scale ensemble
clustering problem is still a significant challenge due to its high
computational complexity, Moreover, it is noteworthy that the ensemble
generation step considerably takes up the run-time during large-scale ensemble
clustering tasks, which has been rarely investigated in the literature.

In light of this, we propose a large-scale ensemble spectral clustering (LSEC)
method to alleviate the problem of the application of ensemble clustering for
large-scale data. In LSEC, a spectral clustering based ensemble generation
method is designed to handle nonlinear datasets efficiently and provide high-
quality base clusterings. The ensemble generation process is further
accelerated by reusing $K$-nearest neighbors among base clusterings and using
light-$k$-means to obtain the clustering results. After ensemble generation, a
bipartite graph between data points and clusters from base clusterings is
constructed to produce consensus clustering through the bipartite graph
partitioning method efficiently. Experimental results on ten large-scale data
sets demonstrate that LSEC delivers highly efficient and high-quality
clustering performance compared to some state-of-the-art consensus clustering
methods.

The main contributions of the proposed method are summaries as follows:

  * •

An efficient spectral clustering based ensemble generation method is designed
to handle large-scale datasets and provide high-quality base clusterings via
divide-and-conquer based large-scale spectral clustering method.

  * •

Two accelerating tricks are proposed: 1) the computation of similarity among
multiple base clusterings is accelerated by reusing the $K$-nearest neighbors;
2) the process of obtaining base clustering results is accelerated by the
light-$k$-means method.

  * •

The proposed method efficiently generates base clusterings and conducts
bipartite graph partitioning to find the consensus clustering. Its
computational and space complexity is dominated by $O(\frac{m}{q}N\alpha d)$
and $O(NK)$, which achieves a lower computational complexity than most
existing ensemble clustering methods.

##  2 Preliminaries

###  2.1 Ensemble Clustering

Ensemble clustering aims to combine multiple base clustering algorithms to
achieve better clustering results. Let ${X}$ be a dataset
${X}=\left\\{x_{1},\ldots,x_{n}\right\\}$ with $n$ data points. The ensemble
generation is the first step, which applies a specific clustering algorithm to
produce $m$ base clusterings. Let
${\Pi}=\left\\{\pi_{1},\ldots,\pi_{m}\right\\}$ be a set of base clusterings,
where $pi_{i}$ is $i$-th base clustering and
$\pi_{i}=\left\\{\pi_{i}\left(x_{1}\right),\pi_{i}\left(x_{2}\right),\cdots,\pi_{i}\left(x_{n}\right)\right\\}$
indicates the clustering labels for all data points. Many studies [6, 14, 32,
19, 26, 9] use $k$-means based ensemble generation while some studies [13, 16]
points out that spectral clustering based ensemble generation can
significantly improve clustering quality on the nonlinear datasets. After
ensemble generation, the consensus function is used to integrate all base
clusterings into a consensus one, which is the second step.

###  2.2 Divide-and-conquer based large-scale spectral clustering algorithm

Divide-and-conquer based large-scale spectral clustering algorithm (DnC-SC)
has been proposed as an effective method for large-scale clustering tasks
[17]. It first constructs an approximate similarity matrix via a divided-and-
conquer based landmark selection and approximates $K$-nearest landmark
searching. Then, it transfers the original spectral clustering problem into a
bipartite graph partition problem to find the low-dimensional embedding by
solving a smaller eigen-problem. Finally, it applies $k$-means on the low-
dimensional embedding to obtain the final clustering result.

Let ${R}=\\{r_{1},r_{2},\cdots,r_{p}\\}$ denote a set of landmarks, where
$r_{i}\in\mathbb{R}^{d}$ has the same dimension as $x_{i}$. The divided-and-
conquer based landmark selection is designed to generate a set of landmark
points which can best represent the original data $X$. The objective function
(1) measure how well $R$ represent $X$ by compute the residual sum of squares
(RRS) between each $x_{j}$ and its nearest $r_{i}$.

$\centering\begin{split}\zeta=\sum_{i=1}^{p}\sum_{x_{j}\in
S_{i}}\left\|x_{j}-r_{i}\right\|^{2},\\\ \end{split}\@add@centering$ (1)

where $\zeta$ denotes RSS and $S_{i},S_{2},\dots,S_{p}$ indicate the subsets
that are nearest to $r_{1},r_{2},\cdots,r_{p}$, respectively. For each subset
$S_{i}$, $r_{i}$ is the subset center. The objective function (1) can be
rewritten as follows:

$g(X,p)=\operatorname*{arg\,min}_{S_{1},\dots,S_{p}}\sum_{i=1}^{p}\sum_{x_{j}\in
S_{i}}\left\|x_{j}-r_{i}\right\|^{2}.$ (2)

The recursive functions (3) and (4) are used to divide the optimization
problem into small sub-problems which are easier to be solved. The parameter
$\alpha$ is used to determine the upper bound of $k_{i}$, which controls the
landmark selection rate.

$\displaystyle g(Q,h)=\bigcup_{i=1}^{m}g(A_{i},k_{i}),$ (3)
$\displaystyle\\{A_{1},\dots,A_{m}\\}=g(Q,m),$ (4)

The light-$k$-means algorithm [17] is used to solve the larger dividing
process $g(\cdot)$ (with more than $10p$ samples), which randomly selects a
part of samples to find subset by $k$-means and then assign remained data
points to the nearest subsets. For the smaller dividing processes (with less
than or equal to $10p$ samples), $k$-means is directly used to find the
subsets.

The similarities between each $x_{i}\in X$ and its $K$-nearest landmarks are
used to construct a sparse similarity matrix. The centers’ nature of landmarks
is used to estimate the $K$-nearest landmarks. Let $S_{x_{i}}$ be the subset
and $x_{i}\in S_{x_{i}}$. Denote $r^{1}_{x_{i}}$ is the landmark that is the
center of $S_{x_{i}}$. According to the center’s nature of landmark,
$r^{1}_{x_{i}}$ is treated as the nearest landmark of $x_{i}$. In DnC-SC, a
set of $K^{\prime}$-nearest landmarks ($K^{\prime}>K$) of $r^{1}_{x_{i}}$ is
first obtained, denoted as ${N}_{K^{\prime}}(r_{x_{i}}^{1})$; then $K$-nearest
landmark of $x_{i}$ are searched from ${N}_{K^{\prime}}(r_{x_{i}}^{1})$,
denoted as ${N}_{K}(x_{i})$. Finally, the sparse similarity matrix $B$ is
constructed as follows [35, 36]:

$\displaystyle b_{ij}$
$\displaystyle=\begin{cases}\exp(\frac{-\left\|x_{i}-r_{j}\right\|^{2}}{2\sigma^{2}}),&\text{if
$r_{j}\in N_{K}(x_{i})$},\\\ 0,&\text{otherwise,}\end{cases}$ (5)

where the Gaussian kernel is used to measure the similarity and $\sigma$ is
the bandwidth parameter.

The similarity matrix $B$ reflects the relationship between data $X$ and
landmarks $R$, which can be treated as the edge of the bipartite graph
$G(X,R,B)$. Therefore, the spectral clustering problem is converted into a
bipartite graph partition problem. According to [17], the low-dimensional
embedding of $R$ side can be computed as follows:

$\displaystyle L_{{R}}V=\lambda D_{{R}}V,$ (6) $\displaystyle U=D_{X}^{-1}BV.$
(7)

where $L_{{R}}=D_{{R}}-B^{T}D_{X}^{-1}B$, $D_{X}\in\mathbb{R}^{n\times n}$ and
$D_{{R}}\in\mathbb{R}^{p\times p}$ are the diagonal matrices whose entries are
$d_{X}(i,i)=\sum_{j=1}^{n}B_{ij}$ and $d_{{R}}(j,j)=\sum_{i=1}^{n}B_{ij}$,
respectively. (6) is a small eigen-problem with size $p\times p$. $U$ is the
$c$ bottom eigenvectors of $X$ side. Finally, $k$-means is conducted on $U$ to
find $c$ clusters as the final clustering result.

##  3 Proposed Framework

To improve the scalability of ensemble clustering, we propose the LSEC method
that complies with the large-scale spectral clustering based formulation and
aims to break through the efficiency bottleneck of previous algorithms. LSEC
method consists of two steps: (1) Large-scale spectral clustering based
ensemble generation: we design a new framework that applies the state-of-the-
art large-scale spectral clustering algorithm to product base clusterings and
further accelerate the process by reusing the $K$-nearest landmarks and using
light-$k$-means to obtain base clustering results. (2) Bipartite graph
partitioning based consensus function: we construct a bipartite graph between
data points and clusters from base clusterings and obtain the consensus
clustering result by bipartite graph partitioning. Fig. 1 shows an overview of
proposed method.

###  3.1 Ensemble Generation based on Large-scale Spectral Clustering

The ensemble generation step aims to produce diverse $m$ base clusterings with
high efficiency. To improve the scalability of ensemble generation, we
consider the divide-and-conquer based large-scale spectral clustering [17] as
the base clustering algorithm, which can better handle nonlinear datasets than
transitional clustering algorithm like $k$-means and maintain high efficiency.
For better diversity of base clusterings and higher efficiency, we construct
similarity matrices with multiple $K$-nearest neighbors graph of
sparsification via reusing the $K$-nearest landmarks. Moreover, the bipartite
graph partitioning is accelerated by applying light-$k$-means to obtain the
clustering results.

Figure 1: An overview of proposed method. Given a dataset, $\frac{m}{q}$ sets
of landmarks are first generated, then a set of $K$-nearest neighbors are
found for each $R^{(i)}$ and $m$ sparse similarity matrices are constructed,
finally the base clusterings are obtained through a bipartite graph
partitioning process. The proposed method accelerates the similarity matrix
construction by recycling $K$-nearest neighbors and bipartite graph
partitioning by applying light-$k$-means.

####  3.1.1 Landmark Selection

First, the $\frac{m}{q}$ sets of landmarks are independently generated by
solving the optimization problem (2). We recursively apply (3) and (4) to find
an approximate local solution and turn the subset centers as landmarks. Let
${R}^{(i)}=\\{r_{1}^{(i)},r_{2}^{(i)},\cdots,r_{p}^{(i)}\\}$ is a set of
landmarks. Repeat the divide-and-conquer based landmark selection
$\frac{m}{q}$ times, we have $\frac{m}{q}$ sets of landmarks as follows:

$\displaystyle\mathcal{R}=\\{{R}^{(1)},{R}^{(2)},\cdots,{R}^{(\frac{m}{q})}\\},$
(8)

where ${R}^{(i)}$ indicates the $i$-th set of landmarks and ${R}$ is a set
containing all ${R}^{(i)}$, $i=1,2,...,\frac{m}{q}$. The generation of each
${R}^{(i)}$ costs $O(N\alpha d)$ time complexity and constructing
$\mathcal{R}$ totally costs $O(\frac{m}{q}N\alpha d)$ time complexity.

####  3.1.2 Searching $K$-nearest landmarks

To construct a sparse similarity matrix with $K$-nearest neighbor
sparsification, we need to search the $K$-nearest landmarks for each data
point $x_{i}$. We consider constructing multiple similarity matrices with
different sparsification of $K$-nearest landmarks for better diversity of base
clusterings. Let $K_{1}<K_{2}<\dots<K_{q}$ be a set of numbers. We search
$K_{1},K_{2},\dots,K_{q}$-nearest landmarks for each $x_{i}$, denoted as
${N}_{K_{1}}(x_{i}),{N}_{K_{2}}(x_{i}),\dots,{N}_{K_{q}}(x_{i})$. According to
the definition of $K$-nearest neighbors, we have

${N}_{K_{1}}(x_{i})\subset{N}_{K_{2}}(x_{i})\subset\dots\subset{N}_{K_{q}}(x_{i}).$
(9)

That is, $N_{K_{j_{1}}}(x_{i})$ is a subset of $N_{K_{j_{2}}}(x_{i})$ if
$K_{j_{1}}<K_{j_{2}}$. Therefore, we only need to compute $K_{q}$-nearest
landmarks and then obtain the other $K_{1},K_{2},\dots,K_{q}$-nearest
landmarks based on it without recomputing. We call this process reusing the
nearest landmarks. Reusing the nearest landmarks accelerates the process of
spectral clustering based ensemble generation. It directly reduces the
computational time in two high-cost steps, landmark selection and searching
$K$-nearest landmarks, by nearly $q$ times. Besides efficiency, it also
enhances the diversity of base clusterings by exploring multiple nearest
neighbor graphs, which is helpful to improve the effectiveness of the proposed
method.

####  3.1.3 Similarity Matrix Construction

Then, the sparse similarity matrix between $X$ and each ${R}^{(i)}$ is
constructed according to (5). Instead of constructing one similarity matrix
for one set of landmarks, we build multiple similarity matrices using
different sets of landmarks. For each ${R}^{(i)}$, we construct $q$ sparse
similarity matrices with $K_{1},K_{2},\dots,K_{q}$-nearest landmarks according
to (5), respectively. We construct $m$ similarity matrices as follows:

$\mathcal{B}=\\{{B}^{(1)}_{1},\dots,{B}^{(1)}_{q},{B}^{(2)}_{1},\dots,{B}^{(2)}_{q},\dots,{B}^{(\frac{m}{q})}_{1},\dots,{B}^{(\frac{m}{q})}_{q}\\},$
(10)

where ${B}^{(i)}_{j}$ indicates a similarity matrix between $X$ and $R^{(i)}$
with sparsification of $K_{j}$-nearest landmarks, $\mathcal{B}$ is a set
containing all ${B}^{(i)}_{j}$ and the total size of $\mathcal{B}$ is $m$. The
computational cost to obtain a sparse similarity matrix ${B}^{(i)_{j}}$ is
$O(NK_{j}d)$ [17]. By reusing the nearest landmarks, we can generate
$\mathcal{B}$ with only $O(\frac{m}{q}NK_{q}d)$ computational cost. For
convenience, we will use $K$ instead of $K_{q}$ to show the computational
complexity in the rest of the paper.

####  3.1.4 Bipartite Graph Partitioning

After obtaining $m$ similarity matrices, we treat each ${B}^{(i)}_{j}$ as the
edge of a bipartite graph $G(X,R^{(\lceil\frac{i}{q}\rceil)},{B}^{(i)_{j}})$
and solve a bipartite graph partition problem by (6) and (7) to construct a
$c^{(i)}_{j}$-dimensional embedding denoted as ${U}^{(i)_{j}}$. Note that
$c^{(i)}_{j}$ is also the number of clusters. It costs $O(p^{3})$ time
complexity to solve each bipartite graph partition problem (6) and
$O(NK(K+c^{(i)}_{j}))$ to compute the each $c^{(i)}_{j}$-dimensional
embedding. The cluster number of $c^{(i)}_{j}$ is randomly selected as follow:

$\displaystyle c^{(i)}_{j}=\lfloor\tau(c_{max}-c_{min})\rfloor+c_{min},$ (11)

where $\tau\in[0,1]$ is a random variable and $c_{max}$ and $c_{min}$ are the
upper and lower bounds of the cluster number, respectively.

The obtained $c^{(i)}_{j}$ eigenvectors are stacked to form a new matrix, upon
which the light-$k$-means [17] is applied to construct the base clustering
result. In light-$k$-means, a set of $p^{\prime}$ samples are first randomly
selected as representatives, then $c$ clusters centers are generated by
applying $k$-means clustering on $p^{\prime}$ representatives, finally, assign
labels to remained samples according to their nearest cluster centers. The
computational complexity of light-$k$-means is $O(pcdt+Ncd)$, where $O(Ncd)$
is the dominated term and $d$ is the dimensional size. The light-$k$-means
alleviates the computational cost from $t$ iterations and can achieve more
efficiency on the platform optimized for matrix operation. The use of
light-$k$-means significantly accelerates the process of obtaining base
clusterings for large-scale datasets. Finally, $m$ base clusterings are
generated, which are represented as

$\displaystyle\Pi=\\{{\pi}^{(1)}_{1},\dots,{\pi}^{(1)}_{q},{\pi}^{(2)}_{1},\dots,{\pi}^{(2)}_{q},\dots,{\pi}^{(\frac{m}{q})}_{1},\dots,{\pi}^{(\frac{m}{q})}_{q}\\},$
(12)

where $\pi^{(i)}_{j}$ denotes a base clustering with $c^{(i)}_{j}$ clusters.
For convenience, we use $c$ instead of $c^{(i)}_{j}$ to show the computational
complexity in the rest paper. The computational complexity of using
light-$k$-means is $O(Nc^{2}+p^{\prime}c^{2}t)$, where $O(Nc^{2})$ is the
dominated term. Overall, the computational complexity of the bipartite graph
partition is $O(m(N(K^{2}+c^{2}+Kc)+p^{3}))$. We summarize the ensemble
generation process of the proposed method in algorithm 1.

1

Input: Dataset $X$, number of base clusterings $m$, a set of number of
$K$-nearest landmarks $K_{1},K_{2},\dots,K_{q}$

Output: base clusterings $\pi_{1},\pi_{2},\dots,\pi_{m}$

2

3Solve (1) by recursively applying (3) and (4) to obtain $\frac{m}{q}$ sets of
landmarks $\mathcal{R}$;

4 for  _$i\leftarrow 1$ to $\frac{m}{q}$_ do

5 Search $K_{q}$-nearest landmarks of each data points according to [17];

6 for  _$j\leftarrow 1$ to $q$_ do

7 Obtain $K_{j}$-nearest landmarks of each data points according to (9);

8 Construct similarity matrix between $X$ and $R^{i}$ with sparsification of
$K_{j}$-nearest landmarks by (5);

9 end for

10

11 end for

12Collect all similarity matrices $\mathcal{B}$ by (10);

13 for  _$i\leftarrow 1$ to $m$_ do

14 Find a low-dimensional embedding $U$ by (6) and (7);

15 Apply light-$k$-means on the embedding $U$ to obtain base clustering
$\pi_{i}$.

16 end for

Algorithm 1 Proposed ensemble generation

###  3.2 Consensus Function based on Bipartite Graph Partitioning

After ensemble generation, the base clusterings will be combined according to
a consensus function for obtaining the consensus partition. Again, we treat
this problem as a bipartite graph partition problem and give a similar
solution like section 2.2.

Table 1: The cluster indicator matrix | $\omega_{1}$ | $\omega_{2}$ | $\cdots$ | $\omega_{C}$ | $\sum$
---|---|---|---|---|---
$\mathcal{X}_{1}$ | $\tilde{b}_{11}$ | $\tilde{b}_{12}$ | $\cdots$ | $\tilde{b}_{1C}$ | $m$
$\mathcal{X}_{2}$ | $\tilde{b}_{21}$ | $\tilde{b}_{22}$ | $\cdots$ | $\tilde{b}_{2C}$ | $m$
$\cdot$ | $\cdot$ | $\cdot$ | $\cdots$ | $\cdot$ | $\cdot$
$\mathcal{X}_{n}$ | $\tilde{b}_{n1}$ | $\tilde{b}_{n2}$ | $\cdots$ | $\tilde{b}_{nC}$ | $m$
$\sum$ | $\|\omega_{1}\|$ | $\|\omega_{2}\|$ | $\cdots$ | $\|\omega_{C}\|$ | $Nm$

To define the bipartite graph, we first collect all clusters though the base
clusterings as (12) and we denotes the clusters in (13) for clarity.

$\displaystyle\Psi=\\{{\Omega}^{(1)}_{1},\dots,{\Omega}^{(1)}_{q},{\Omega}^{(2)}_{1},\dots,{\Omega}^{(2)}_{q},\dots,{\Omega}^{(\frac{m}{q})}_{1},\dots,{\Omega}^{(\frac{m}{q})}_{q}\\},$
(13)

where $\Omega^{(i)}_{j}$ indicate the set of clusters in $\pi^{(i)}_{j}$.
There are $c^{(i)}_{j}$ clusters in each $\Omega^{(i)}_{j}$, which we denote
as:

$\Omega^{(i)}_{j}=\\{\omega^{\prime}_{1},\omega^{\prime}_{2},\dots,\omega^{\prime}_{c^{(i)}_{j}}\\},$
(14)

where $\omega^{\prime}_{t}$ is the $t$-th cluster in $\Omega^{(i)}_{j}$. Thus,
the total number of clusters in $\Psi$ can be counted as
$C=\sum_{i=1}^{\frac{m}{q}}\sum_{j=1}^{q}c_{j}^{(i)}$. For convenience, we
simplify the notation of (15) as follows:

$\Psi=\\{\omega_{1},\omega_{2},\dots,\omega_{C}\\},$ (15)

After the definition of $\Omega$, we design a bipartite graph between data
points and clusters as follow:

$\displaystyle\tilde{G}=\\{\mathcal{X},\Omega,\tilde{B}\\},$ (16)

where $\tilde{B}$ is the cross-affinity matrix between $\mathcal{X}$ and
$\Omega$. $\tilde{B}$ can also be interpreted as the cluster indicator matrix
of ${X}$. Table 1 shows the cluster indicator matrix, where $b_{ij}=1$
indicates that ${X}_{i}\in\omega_{j}$. $\tilde{G}$ is an unweighted bipartite
graph where any edge between node $X_{i}$ and $\omega_{j}$ indicates the
cluster relationship $X_{i}\in\omega_{j}$. We can give the formula of
$\tilde{B}$ as follow:

$\displaystyle\tilde{b}_{ij}$ $\displaystyle=\begin{cases}1,&\text{if
$x_{i}\in\omega_{j}$},\\\ 0,&\text{otherwise.}\end{cases}$ (17)

As Table 1 shows, the sum of each rows of $\tilde{B}$ is as the same as number
of base clusterings $m$ because there is only one cluster $x_{i}$ belongs to
each base clustering $\pi^{j}$, i.e., $\forall i^{\prime}\neq j^{\prime}$, if
$\omega_{i^{\prime}}\in\pi^{i}$ and $\omega_{j^{\prime}}\in\pi^{i}$, then
$\omega_{i^{\prime}}\bigcap\omega_{j^{\prime}}=\emptyset$. Though the number
of samples in each $\omega_{i}$ is uncertain, i.e., $\|\omega_{i}\|$, the
total number of non-zeros entries is clearly $Nm$ (see Table 1).

For this modified bipartite graph $\tilde{G}$, we consider a similar partition
strategy to what we introduced in Section 2.2. According to [17], we can write
the full similarity of $\mathcal{G}$ as follow

$\tilde{W}=\left[\begin{array}[]{ll}{0}&\tilde{B}\\\
\tilde{B}^{T}&{0}\end{array}\right].$ (18)

Then the we have a generalized eigen-problem of $\tilde{G}$

$\displaystyle\tilde{L}\tilde{f}=\lambda\tilde{D}\tilde{f},$ (19)

where $\tilde{L}=\tilde{D}-\tilde{W}$ and $\tilde{D}$ is a diagonal matrix
with $\tilde{d_{ii}}=\sum_{j=1}^{n}\tilde{w_{ij}}$. According to (20) and (7),
we design a smaller eigen-problem to compute the eigenvetor $\tilde{U}$ in
$\mathcal{X}$ side as follows:

$\displaystyle{L}_{\Omega}\tilde{V}=\tilde{\lambda}{D}_{\Omega}\tilde{V},$
(20)

where
$L_{\Omega}=\tilde{D}_{\Omega}-\tilde{B}^{\top}{\tilde{D}_{\mathcal{X}}}^{-1}\tilde{B}$
is the graph Laplacian, $\tilde{D}_{\mathcal{X}}\in R^{n\times n}$ and
$\tilde{D}_{\mathcal{R}}\in R^{p\times p}$ are the diagonal matrices whose
entries are $\tilde{d}_{\mathcal{X}}(i,i)=\sum_{j=1}^{n}\tilde{B}_{ij}$ and
$\tilde{d}_{\mathcal{R}}(j,j)=\sum_{i=1}^{n}\tilde{B}_{ij}$, respectively. The
size of $L_{\Omega}$ is $C\times C$. Solving the eigen-problem (20) cost
$O(C^{3})$ computational time. Substituting $\tilde{V}$ into (7), we can
computer $\tilde{U}$ as follow

$\tilde{D}=\tilde{D}_{\mathcal{X}}^{-1}\tilde{B}\tilde{V}.$ (21)

The $\tilde{c}$ bottom eigenvectors $\tilde{U}$ can be computed with with
$O(Nm(m+c))$ time. Finally, the consensus clustering results in LSEC can be
obtained by the $k$-means method with $O(Nc^{2}t)$ time. We summarize the
proposed method LSEC in algorithm 2.

Input: Dataset $X$, number of base clusterings $m$, a set of number of
$K$-nearest landmarks $K_{1},K_{2},\dots,K_{q}$, $m$ base clusterings, number
of clusters $\tilde{c}$

Output: Consensus clustering $\tilde{\pi}$

1 Produce $m$ base clustering by large-scale ensemble generation;

2 Construct the cluster indicator matrix $\tilde{B}$ according to (17);

3 Solve the eigen-problem (20) to compute $\tilde{V}$;

4 Find a low-dimensional embedding $\tilde{U}$ of $X$ by (21);

Applying $k$-means to find $\tilde{c}$ clusters on $\tilde{U}$ to obtain
consensus clustering. Obtain consensus clustering by large-scale consensus
function.

Algorithm 2 Large-scale ensemble spectral clustering Table 2: Comparison of
the computational complexity between LSEC and U-SPEC.

Method | Ensemble Generation | Consensus Function
---|---|---
Landmark selection | Similarity construction | Bipartite graph partitioning
U-SPEC | $O(mp^{2}dt)$ | $O(mNp^{\frac{1}{2}}d)$ | $O(m(N(K^{2}+c^{2}t+Kc)+p^{3}))$ | $O(N(m^{2}+mk+c^{2}t)+C^{3})$
LSEC | $O(\frac{m}{q}N\alpha d)$ | $O(\frac{m}{q}NKd)$ | $O(m(N(K^{2}+c^{2}+Kc)+p^{3}))$ | $O(N(m^{2}+mk+c^{2}t)+C^{3})$

##  4 Discussion

###  4.1 Computational Complexity Analysis

In this section, we summarize the computational cost of the proposed method.
The ensemble generation of LSEC algorithm takes $O(mN(\alpha
d+K^{2}+Kc+Kd+qc^{2})+p^{3}+p^{2}(d+K))$ computational cost. The consensus
function of LSEC takes $O(N((qm)^{2}+qmk+c^{2}t)+{C}^{3})$ time. With
consideration to $m,q,k,K<\alpha\ll p\ll N$, the dominant term of the overall
time complexity of LSEC is $O(Nm(\alpha d+qk^{2}))$. Meanwhile, the memory
costs of the ensemble generation and the consensus function of our LSEC
algorithm are respectively $O(N\alpha)$ and $O(Nm)$. Table 2 provides a
comparison of the computational complexity of our DnC-SC algorithm against a
state-of-the-art large-scale ensemble clustering method U-SPEC.

###  4.2 Relations with Other Methods

As a large-scale spectral ensemble clustering method, the proposed method is
closely related to the U-SENC method in [13]. We compare the proposed method
with the U-SENC to discuss the improvements of the proposed method.

Firstly, we compare them to the ensemble generation methods in the term of
diversity and efficiency. In U-SENC, base clusterings are directly generated
using a large-scale spectral clustering U-SPEC with different numbers of
clusters. As a large-scale spectral clustering method, the U-SPEC method also
uses the landmark selection technique. Thus, the diversity of base clusterings
of U-SENC is from two facts: the different landmarks and the number of
clusters of ensemble generation. However, the $K$-nearest neighbor graph is
not used to improve the diversity in U-SENC further. In our proposed method,
we consider the various landmarks and number of clusters and use different
$K$-nearest neighbors to construct a sparse similarity matrix to improve the
overall diversity of base clustering.

Secondly, we compare them to the ensemble generation methods in the term of
efficiency. Since the different $K$-nearest neighbor graphs can share the same
$K$-nearest neighbors between data points and landmarks, the computational
complexity of similarity matrix construction is much less than the U-SENC
method. For large-scale datasets, another computational bottleneck is the
final $k$-means step of large-scale spectral clustering. In our proposed
method, we use the light-$k$-means to accelerate the base clustering results,
significantly improving large-scale datasets’ efficiency.

Overall, LSEC redesigns the ensemble generation framework based on a more
efficient clustering method (i.e., DnC-SC) and accelerates the process by
reusing the $K$-nearest neighbors among multiple base clusterings.
Furthermore, a light-$k$-means method is used to fast obtain the base
clustering results. The computational complexity of the proposed method is
faster than most existing large-scale ensemble clustering methods.

##  5 Experiments

In this section, we conduct experiments on five real and five synthetic
datasets to evaluate the performance of the proposed LSEC method. The
comparison experiments against several state-of-the-art spectral clustering
methods show better clustering quality and efficiency for LSEC methods.
Besides that, the analysis of parameters is performed. For each experiment,
the test method is repeated 20 times, and the average performance is reported.
All experiments are conducted in Matlab R2020a on a Mac Pro with 3 GHz 8-Core
Intel Xeon E5 and 16 GB of RAM.

###  5.1 Datasets and Evaluation Measures

Table 3: Properties of the real and synthetic datasets. Dataset |  #Object |  #Dimension |  #Class
---|---|---|---
|  _USPS_ |  9298 |  256 |  10
_Real_ |  _PenDigits_ |  10,992 |  16 |  10
|  _Letters_ |  20,000 |  16 |  26
|  _MNIST_ |  70,000 |  784 |  10
|  _Covertype_ |  581,012 |  54 |  7
_Synthetic_ |  _TB-1M_ |  1,000,000 |  2 |  3
|  _SF-2M_ |  2,000,000 |  2 |  4
|  _CC-5M_ |  5,000,000 |  2 |  3
|  _CG-10M_ |  10,000,000 |  2 |  11
|  _FL-20M_ |  20,000,000 |  2 |  13

(a) _TB-1M_ ($0.1\%$) (b) _SF-2M_ ($0.1\%$)

(c) _CC-5M_ ($0.1\%$) (d) _CG-10M_ ($0.1\%$)

(e) _FL-20M_ ($0.1\%$)

Figure 2: Illustration of the five synthetic datasets. Note that only $0.1\%$
samples of each dataset are plotted.

Our experiments are conducted on ten large-scale datasets, varying from nine
thousands to as large as twenty million data points. Specifically, the five
real datasets are _PenDigits_ [1]
111https://archive.ics.uci.edu/ml/datasets/Pen-
Based+Recognition+of+Handwritten+Digits, _USPS_ [4] 222
http://www.cad.zju.edu.cn/home/dengcai/Data/MLData.html, _Letters_ [7]
333https://archive.ics.uci.edu/ml/datasets/Letter+Recognition, _MNIST_ [3],
and _Covertype_ [2] 444https://archive.ics.uci.edu/ml/datasets/covertype. The
five synthetic datasets are _Two Bananas_ (_TB-1M_), _Smiling Face-2M_
(_SF-2M_), _Concentric Circles-5M_ (_CC-5M_), _Circles and Gaussians-10M_
(_CG-10M_), _Flower-20M_ (_FL-20M_) [13]
555https://www.researchgate.net/publication/330760669. Fig. 2 shows the
synthetic datasets. The properties of the datasets are summarized in Table 3.

We adopt two widely used evaluation metrics, i.e., Normalized Mutual
Information (NMI) [23] and Accuracy (ACC) [34], to evaluate the clustering
results. Let $X=[x_{1},x_{2},...,x_{n}]$ be the data matrix. For each data
point $x_{i}$, denote $\pi_{t}(x_{i})$ and $\pi_{c}(x_{i})$ as the cluster
label of ground truth and obtained cluster label from clustering methods,
respectively. The ACC is defined as:

$\text{ACC}=\frac{\sum_{i=1}^{n}\delta(\pi_{t}(x_{i}),\text{map}(\pi_{c}(x_{i})))}{n},$
(22)

where $n$ is the number of data and $\delta(\pi_{t}(x_{i}),\pi_{c}(x_{i}))$ is
a function to check $\pi_{t}(x_{i})$ and $\pi_{c}(x_{i})$ are equal or not,
returning 1 if equals otherwise returning 0. The map$(\pi_{c}(x_{i}))$ is a
best mapping function that maps each predicted label to the most possibly true
cluster label by permuting operations [33].

The NMI is the normalization of Mutual information by the joint entropy as
follow:

$NMI=\frac{\sum_{\pi_{t}(x_{i})\in T,\pi_{c}(x_{i})\in
C}p(\pi_{t}(x_{i}),\pi_{c}(x_{i}))\text{ln}\frac{p(\pi_{t}(x_{i}),\pi_{c}(x_{i}))}{p(\pi_{t}(x_{i}))p(\pi_{c}(x_{i}))}}{-\sum_{\pi_{t}(x_{i})\in
T,\pi_{c}(x_{i})\in
C}p(\pi_{t}(x_{i}),\pi_{c}(x_{i}))\text{ln}(p(\pi_{t}(x_{i}),\pi_{c}(x_{i})))},$
(23)

A better clustering result will provide a larger value of NMI/ACC. Both NMI
and ACC are in the range of $[0,1]$.

###  5.2 Compared Methods and Experimental Settings

In this experiments, we compare the proposed method with one baseline
clustering method, i.e., the divide-and-conquer based large-scale spectral
clustering (DnC-SC) [17], as well as seven state-of-the-art large-scale
spectral clustering methods. The compared spectral clustering methods are
listed as follows:

  1. 1.

EAC [6]: evidence accumulation clustering.

  2. 2.

KCC [32]: $k$-means based consensus clustering.

  3. 3.

PTGP [9]: probability trajectory based graph partitioning.

  4. 4.

SEC [19]: spectral ensemble clustering.

  5. 5.

LWEA [11]: locally weighted evidence accumulation.

  6. 6.

LWGP [11]: locally weighted graph partitioning.

  7. 7.

U-SENC [13]: ultra-scalable ensemble clustering.

Table 4: ACC(%) scores (over 20 runs) by our methods and the baseline ensemble
clustering methods (The best score in each row is in bold).

Dataset | DnC-SC | EAC | KCC | PTGP | SEC | LWGP | U-SENC | LSEC
---|---|---|---|---|---|---|---|---
_PenDigits_ | 82.27±1.33 | 77.67±2.30 | 44.68±5.10 | 80.89±1.26 | 32.37±3.88 | 73.66±2.14 | 87.02±1.65 | 88.26±2.10
_USPS_ | 82.55±1.96 | 66.76±1.69 | 56.37±3.41 | 67.14±0.26 | 40.77±5.82 | 65.82±3.41 | 78.25±2.39 | 80.97±5.31
_Letters_ | 33.54±1.21 | 29.79±0.60 | 24.46±1.24 | 26.66±1.42 | 23.60±1.28 | 27.88±0.78 | 37.03±1.28 | 36.33±0.89
_MINST_ | 74.24±2.14 | N/A | 45.61±4.96 | 66.96±0.68 | 33.15±2.07 | 56.27±1.47 | 75.48±3.01 | 80.19±3.74
_Covertype_ | 23.48±1.86 | N/A | 32.52±0.41 | 23.45±0.96 | 39.63±6.15 | 30.64±0.42 | 21.34±1.06 | 23.42±1.86
_TB-1M_ | 99.62±0.02 | N/A | 67.76±1.41 | 81.95±0.00 | 67.94±3.66 | 99.71±0.45 | 99.75±0.01 | 99.72±2.31
_SF-2M_ | 9.43±0.31 | N/A | 50.94±4.15 | 60.25±0.94 | 49.88±5.68 | 80.04±3.45 | 76.54±2.88 | 85.17±9.18
_CC-5M_ | 99.98±0.00 | N/A | 72.25±6.41 | 34.95±0.00 | 41.57±0.81 | 97.84±3.71 | 99.99±0.00 | 95.61±11.67
_CG-10M_ | 66.83±4.46 | N/A | 58.12±5.41 | 60.89±1.73 | 46.26±5.72 | 71.95±3.19 | 82.34±5.59 | 97.57±3.49
_FL-20M_ | 81.90±5.61 | N/A | 48.21±4.14 | 51.21±1.41 | 41.70±0.42 | 72.15±2.45 | 78.16±3.21 | 82.81±3.21
Avg. score | - | N/A | 50.09 | 55.44 | 41.69 | 67.60 | 72.59 | 77.01
Avg. rank | - | 6.00 | 5.00 | 4.00 | 5.60 | 3.30 | 2.30 | 1.90

Table 5: NMI(%) scores (over 20 runs) by our methods and the baseline ensemble
clustering methods (The best score in each row is in bold).

Dataset | DnC-SC | EAC | KCC | PTGP | SEC | LWGP | U-SENC | LSEC
---|---|---|---|---|---|---|---|---
_PenDigits_ | 82.01±0.21 | 76.12±0.00 | 53.52±3.46 | 78.31±0.34 | 46.44±2.10 | 76.46±1.43 | 83.24±1.11 | 84.65±1.81
_USPS_ | 82.86±1.08 | 69.05±0.00 | 58.27±0.25 | 70.32±1.02 | 49.68±1.89 | 70.71±1.46 | 82.09±1.65 | 83.51±1.34
_Letters_ | 45.37±0.85 | 39.28±0.00 | 34.61±1.40 | 36.98±0.99 | 32.30±0.89 | 39.29±0.46 | 46.40±0.20 | 48.71±0.63
_MINST_ | 72.00±0.51 | N/A | 46.43±4.85 | 62.22±1.12 | 38.84±1.44 | 62.34±0.62 | 75.11±0.58 | 79.42±1.45
_Covertype_ | 8.30±0.30 | N/A | 6.38±3.41 | 8.25±0.43 | 9.23±6.48 | 9.06±0.41 | 9.34±1.21 | 11.64±1.76
_TB-1M_ | 96.42±0.18 | N/A | 24.54±2.45 | 31.89±0.00 | 24.74±4.45 | 97.16±2.41 | 97.45±0.04 | 97.19±9.43
_SF-2M_ | 81.24±0.32 | N/A | 38.06±2.45 | 49.74±0.18 | 33.65±3.22 | 81.95±4.15 | 77.57±2.12 | 84.88±6.55
_CC-5M_ | 99.78±0.01 | N/A | 59.24±0.41 | 0.13±0.00 | 12.93±1.80 | 98.15±7.41 | 99.91±0.00 | 95.57±10.70
_CG-10M_ | 80.91±3.59 | N/A | 63.56±0.41 | 65.09±0.92 | 55.77±6.84 | 78.41±2.93 | 86.28±2.30 | 95.25±1.32
_FL-20M_ | 87.67±3.18 | N/A | 68.10±2.41 | 71.32±1.29 | 53.77±2.52 | 78.51±1.97 | 90.38±2.45 | 91.32±2.44
Avg. score | - | N/A | 45.27 | 47.43 | 35.74 | 69.20 | 74.78 | 77.21
Avg. rank | - | 6.30 | 5.40 | 4.30 | 5.80 | 3.00 | 1.90 | 1.40

Table 6: Time costs(s) of our methods and the baseline ensemble clustering
methods.

Dataset | DnC-SC | EAC | KCC | PTGP | SEC | LWGP | U-SENC | LSEC
---|---|---|---|---|---|---|---|---
_PenDigits_ | 0.64 | 18.78 | 9.19 | 6.06 | 3.00 | 4.00 | 18.31 | 3.40
_USPS_ | 1.25 | 25.79 | 23.56 | 41.32 | 15.08 | 15.93 | 28.82 | 5.81
_Letters_ | 0.90 | 115 | 48.48 | 89.88 | 10.76 | 11.15 | 20.86 | 3.71
_MINST_ | 5.11 | N/A | 831.12 | 2297.05 | 730.33 | 731.64 | 103.35 | 21.36
_Covertype_ | 13.15 | N/A | 634.18 | 16271.2 | 714.86 | 730.33 | 143.41 | 40.16
_TB-1M_ | 5.06 | N/A | 984.15 | 849.62 | 693.67 | 709.60 | 265.80 | 62.50
_SF-2M_ | 13.77 | N/A | 2225.64 | 1475.08 | 1344.66 | 1566.8 | 623.26 | 131.67
_CC-5M_ | 25.37 | N/A | 8541.13 | 3040.33 | 3232.06 | 3006 | 1851.2 | 321.71
_CG-10M_ | 281.05 | N/A | 12351.2 | 7244.01 | 7607.84 | 6685.8 | 3561.4 | 769.51
_FL-20M_ | 837.38 | N/A | 17112.1 | 13343.3 | 14938.73 | 13091 | 11763.07 | 2396.85
Avg. score | - | N/A | 4276.07 | 4465.78 | 2929.10 | 2655.23 | 1837.95 | 375.65
Avg. rank | - | 6.80 | 5.20 | 5.00 | 3.30 | 3.60 | 3.00 | 1.10

There are several common parameters among the methods mentioned above. We set
these parameters as follow:

  * •

We set the number of landmarks as $p=1000$ for LSEC and U-SPEC. The parameter
analysis on $p$ has been conducted in [17].

  * •

We set the $K=5$ for the number of nearest neighbors for LSEC and U-SPEC.

  * •

The DnC-SC method has a unique parameter $\alpha$. In the experiments,
$\alpha=50$ is used for all datasets.

  * •

The base clusterings are generated by $k$-means or large-scale spectral
clustering as suggested by their papers [6, 32, 9, 11, 13]. The number of
cluster $c$ of base clusterings is randomly selected from $[20,60]$. The
number of base clusterings is set as $m=20$ for comparison. The parameter
analysis on $m$ will be shown in Section 5.4.

  * •

The true number of classes on each dataset is used to conduct all experiments.

  * •

Other parameters in the baseline methods are set as suggested by the original
papers.

###  5.3 Comparison Results

Table 7: Clustering performance (ACC(%), NMI(%), and time costs(s)) for different methods by varying number of base clusterings $m$. _Dataset_ |  _MNIST_ |  _Covertype_ |  _TB-1M_ |  _SF-2M_
---|---|---|---|---
ACC |  |  |  |
NMI |  |  |  |
Time cost |  |  |  |
|
Table 8: Clustering performance (ACC(%), NMI(%), and time costs(s)) for LSEC with or without reusing of nearest landmarks. _Data_ |  _MNIST_ |  _Covertype_ |  _TB-1M_ |  _SF-2M_
---|---|---|---|---
ACC |  |  |  |
NMI |  |  |  |
Time cost |  |  |  |
|
Table 9: Clustering performance (ACC(%), NMI(%), and time costs(s)) for LSEC using light-$k$-means or using $k$-means to obtain base clusterings in the ensemble generation. _Data_ |  _MNIST_ |  _Covertype_ |  _TB-1M_ |  _SF-2M_
---|---|---|---|---
ACC |  |  |  |
NMI |  |  |  |
Time cost |  |  |  |
|

The experimental comparison results are reported in Tables 4, 5, and 6. Note
that DnC-SC is not an ensemble clustering algorithm; its clustering results
are provided for reference only.

As shown in Tables 4 and 5, our LSEC algorithm obtains the highest ACC and NMI
scores on most of datasets. In terms of average score across the ten datasets,
LSEC achieves the best average ACC($\%$) and NMI($\%$) scores of $77.01$ and
$77.21$, respectively. While the second-best ensemble clustering method (i.e.,
U-SENC) achieves average ACC($\%$) and NMI($\%$) scores of $72.59$ and
$74.78$, respectively. The EAC, KCC, PTGP, SEC, LWGP methods use the $k$-means
based ensemble generation method. The LSEC and U-SENC methods that use the
spectral clustering based ensemble generation show better clustering quality
of ACC and NMI than others on most datasets. In terms of average rank, LSEC
obtains an average rank of 1.90 w.r.t ACC and 1.40 w.r.t. NMI, while the
second-best method obtains an average rank of 2.30 w.r.t. ACC and 1.90 w.r.t.
NMI.

In Table 6, the time costs of different ensemble clustering methods are
provided. The proposed LSEC method achieves the lowest time costs on nine
datasets and the second-lowest time cost on one dataset. Except _PenDigits_
dataset, LSEC is 2.4 (_FL-20M_) to 5.75 (_CC-5M_) times ahead of the second-
best method in time consumption. The LSEC method has shown its significant
advantage over other ensemble clustering methods, especially on large-scale
datasets.

###  5.4 Parameter analysis on Ensemble Size $m$

We conduct a parameter analysis experiment to demonstrate the performance of
the proposed method, varying different parameter values of $m$. The parameter
$m$ denotes the number of base clusterings, which is a common parameter in all
ensemble clustering methods. We select four dataset (_MNIST_ , _Covertype_ ,
_TB-1M_ and _SF-2M_) as benchmark datasets to conduct the following
experiments. As shown in Table 7, LSEC shows better performance of ACC and NMI
than most other ensemble clustering methods except ACC score on _Covertype_
dataset. Meanwhile, LSEC consistently requires a lower computational cost than
all other ensemble clustering methods.

###  5.5 Influence of reusing of $K$-nearest landmarks

In this section, we compare the performances of the proposed method with or
without reusing of nearest landmarks, denoted as LSEC and LSEC-without-
reusing. The experimental results are reported in Table 8. As we mentioned,
the reusing of nearest landmarks brings better efficiency in searching
$K$-nearest landmarks. In Table 8, LSEC and LSEC-without-reusing show similar
performances to each other, but LSEC cost obviously less time. Since reusing
of nearest landmarks does not influence the accuracy of nearest landmarks, we
consider that the difference of ACC and NMI between the two methods comes from
the randomness of the algorithm. This result indicates that reusing of nearest
landmarks achieves significantly better efficiency while maintaining a similar
clustering result.

###  5.6 Influence of light-$k$-means

In this section, we compare the performances of the proposed method using
light-$k$-means or using $k$-means to obtain base clusterings in the ensemble
generation. The experimental results are reported in Table 9. Generally, two
methods show the similar performance of ACC and NMI. Especially, LSEC achieves
slightly better ACC and NMI on _MNIST_ and _SF-2F_ datasets, which is possible
because the light-$k$-means method can provide better diversity of base
clusterings on these datasets. Overall, using light-$k$-means in ensemble
generation significantly improves the efficiency of LSEC and yields similar
clustering quality compared to the $k$-means method.

##  6 Conclusion

In this paper, we propose a large-scale spectral ensemble clustering (LSEC)
method to strike a better balance between the efficiency and effectiveness of
ensemble clustering on large-scale datasets. We design an efficient ensemble
generation framework to produce based clustering, applying divide-and-conquer
large-scale spectral clustering to find high-quality base clusterings. In the
ensemble generation of the proposed method, we accelerate the process of
searching $K$-nearest neighbors by reusing strategy and obtaining base
clustering by the light-$k$-means method. After the ensemble generation step,
we combine all based clustering into a consensus cluster through a bipartite
graph partitioning based consensus function. The proposed method achieves
lower computational complexity than most existing ensemble clustering methods.
Experiments conducted on ten large-scale datasets show that the proposed
method outperforms other state-of-the-art large-scale spectral clustering
methods.

## Acknowledgment

This study was supported by the New Energy and Industrial Technology
Development Organization (NEDO) Grant (ID:18065620) and JST COI-NEXT.

## References

  * [1] Arthur Asuncion and David Newman.  Uci machine learning repository, 2007.
  * [2] Jock A Blackard and Denis J Dean.  Comparative accuracies of artificial neural networks and discriminant analysis in predicting forest cover types from cartographic variables.  Computers and electronics in agriculture, 24(3):131–151, 1999.
  * [3] Deng Cai, Xiaofei He, and Jiawei Han.  Speed up kernel discriminant analysis.  The VLDB Journal, 20(1):21–33, 2011.
  * [4] Deng Cai, Xiaofei He, Jiawei Han, and Thomas S Huang.  Graph regularized nonnegative matrix factorization for data representation.  IEEE transactions on pattern analysis and machine intelligence, 33(8):1548–1560, 2010.
  * [5] Xiaoli Zhang Fern and Carla E Brodley.  Solving cluster ensemble problems by bipartite graph partitioning.  In Proceedings of the twenty-first international conference on Machine learning, page 36, 2004.
  * [6] Ana LN Fred and Anil K Jain.  Combining multiple clusterings using evidence accumulation.  IEEE transactions on pattern analysis and machine intelligence, 27(6):835–850, 2005.
  * [7] Peter W Frey and David J Slate.  Letter recognition using holland-style adaptive classifiers.  Machine learning, 6(2):161–182, 1991.
  * [8] Dong Huang, Jian-Huang Lai, and Chang-Dong Wang.  Combining multiple clusterings via crowd agreement estimation and multi-granularity link analysis.  Neurocomputing, 170:240–250, 2015.
  * [9] Dong Huang, Jian-Huang Lai, and Chang-Dong Wang.  Robust ensemble clustering using probability trajectories.  IEEE transactions on knowledge and data engineering, 28(5):1312–1326, 2015.
  * [10] Dong Huang, Jianhuang Lai, and Chang-Dong Wang.  Ensemble clustering using factor graph.  Pattern Recognition, 50:131–142, 2016.
  * [11] Dong Huang, Chang-Dong Wang, and Jian-Huang Lai.  Locally weighted ensemble clustering.  IEEE transactions on cybernetics, 48(5):1460–1473, 2017.
  * [12] Dong Huang, Chang-Dong Wang, Hongxing Peng, Jianhuang Lai, and Chee-Keong Kwoh.  Enhanced ensemble clustering via fast propagation of cluster-wise similarities.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2018\.
  * [13] Dong Huang, Chang-Dong Wang, Jian-Sheng Wu, Jian-Huang Lai, and Chee-Keong Kwoh.  Ultra-scalable spectral clustering and ensemble clustering.  IEEE Transactions on Knowledge and Data Engineering, 32(6):1212–1226, 2019.
  * [14] Natthakan Iam-On, Tossapon Boongeon, Simon Garrett, and Chris Price.  A link-based cluster ensemble approach for categorical data clustering.  IEEE Transactions on knowledge and data engineering, 24(3):413–425, 2010.
  * [15] Vladimir Yu Kiselev, Kristina Kirschner, Michael T Schaub, Tallulah Andrews, Andrew Yiu, Tamir Chandra, Kedar N Natarajan, Wolf Reik, Mauricio Barahona, Anthony R Green, et al.  Sc3: consensus clustering of single-cell rna-seq data.  Nature methods, 14(5):483, 2017.
  * [16] Hongmin Li, Xiucai Ye, Akira Imakura, and Tetsuya Sakurai.  Ensemble learning for spectral clustering.  In 2020 IEEE International Conference on Data Mining (ICDM), pages 1094–1099. IEEE, 2020.
  * [17] Hongmin Li, Xiucai Ye, Akira Imakura, and Tetsuya Sakurai.  Divide-and-conquer based large-scale spectral clustering.  arXiv preprint, page 2104.15042, 2021.
  * [18] Tao Li, Mitsunori Ogihara, and Sheng Ma.  On combining multiple clusterings.  In Proceedings of the thirteenth ACM international conference on Information and knowledge management, pages 294–303, 2004.
  * [19] Hongfu Liu, Tongliang Liu, Junjie Wu, Dacheng Tao, and Yun Fu.  Spectral ensemble clustering.  In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pages 715–724, 2015.
  * [20] Hongfu Liu, Junjie Wu, Tongliang Liu, Dacheng Tao, and Yun Fu.  Spectral ensemble clustering via weighted k-means: Theoretical and practical evidence.  IEEE transactions on knowledge and data engineering, 29(5):1129–1143, 2017.
  * [21] Hongfu Liu, Rui Zhao, Hongsheng Fang, Feixiong Cheng, Yun Fu, and Yang-Yu Liu.  Entropy-based consensus clustering for patient stratification.  Bioinformatics, 33(17):2691–2698, 2017.
  * [22] Murilo Coelho Naldi, ACPLF Carvalho, and Ricardo JGB Campello.  Cluster ensemble selection based on relative validity indexes.  Data Mining and Knowledge Discovery, 27(2):259–289, 2013.
  * [23] Noam Slonim and Naftali Tishby.  Agglomerative information bottleneck.  In Advances in neural information processing systems, pages 617–623, 2000.
  * [24] Alexander Strehl and Joydeep Ghosh.  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  Journal of machine learning research, 3(Dec):583–617, 2002.
  * [25] Aditya Tandon, Aiiad Albeshri, Vijey Thayananthan, Wadee Alhalabi, and Santo Fortunato.  Fast consensus clustering in complex networks.  Physical Review E, 99(4):042301, 2019.
  * [26] Zhiqiang Tao, Hongfu Liu, Sheng Li, and Yun Fu.  Robust spectral ensemble clustering.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, pages 367–376. ACM, 2016.
  * [27] Alexander Topchy, Anil K Jain, and William Punch.  Combining multiple weak clusterings.  In Third IEEE international conference on data mining, pages 331–338. IEEE, 2003.
  * [28] Sandro Vega-Pons and José Ruiz-Shulcloper.  A survey of clustering ensemble algorithms.  International Journal of Pattern Recognition and Artificial Intelligence, 25(03):337–372, 2011.
  * [29] Chen Wang, Sebastian M Armasu, Kimberly R Kalli, Matthew J Maurer, Ethan P Heinzen, Gary L Keeney, William A Cliby, Ann L Oberg, Scott H Kaufmann, and Ellen L Goode.  Pooled clustering of high-grade serous ovarian cancer gene expression leads to novel consensus subtypes associated with survival and surgical outcomes.  Clinical Cancer Research, 23(15):4077–4085, 2017.
  * [30] Fei Wang, Xin Wang, and Tao Li.  Generalized cluster aggregation.  In Twenty-First International Joint Conference on Artificial Intelligence, 2009.
  * [31] Xi Wang, Chunyu Yang, and Jie Zhou.  Clustering aggregation by probability accumulation.  Pattern Recognition, 42(5):668–675, 2009.
  * [32] Junjie Wu, Hongfu Liu, Hui Xiong, Jie Cao, and Jian Chen.  K-means-based consensus clustering: A unified view.  IEEE transactions on knowledge and data engineering, 27(1):155–169, 2014.
  * [33] Wei Xu, Xin Liu, and Yihong Gong.  Document clustering based on non-negative matrix factorization.  In Proceedings of the 26th annual international ACM SIGIR conference on Research and development in informaion retrieval, pages 267–273, 2003.
  * [34] Donghui Yan, Ling Huang, and Michael I Jordan.  Fast approximate spectral clustering.  In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 907–916. ACM, 2009.
  * [35] Xiucai Ye and Tetsuya Sakurai.  Robust similarity measure for spectral clustering based on shared neighbors.  ETRI journal, 38(3):540–550, 2016.
  * [36] Xiucai Ye and Tetsuya Sakurai.  Spectral clustering with adaptive similarity measure in kernel space.  Intelligent Data Analysis, 22(4):751–765, 2018.
  * [37] Li Zheng, Tao Li, and Chris Ding.  A framework for hierarchical ensemble clustering.  ACM Transactions on Knowledge Discovery from Data (TKDD), 9(2):1–23, 2014.
  * [38] Caiming Zhong, Xiaodong Yue, Zehua Zhang, and Jingsheng Lei.  A clustering ensemble: Two-level-refined co-association matrix with path-based transformation.  Pattern Recognition, 48(8):2699–2709, 2015.