IEEE TRANSACTIONS ON XXXX

1

Non-iterative Label Propagation on
Optimal Leading Forest

arXiv:1709.08426v1 [cs.LG] 25 Sep 2017

Ji Xu, Guoyin Wang, Senior Member, IEEE

Abstract—Graph based semi-supervised learning (GSSL) has
intuitive representation and can be improved by exploiting
the matrix calculation. However, it has to perform iterative
optimization to achieve a preset objective, which usually leads to
low efficiency. Another inconvenience lying in GSSL is that when
new data come, the graph construction and the optimization have
to be conducted all over again. We propose a sound assumption,
arguing that: the neighboring data points are not in peer-to-peer
relation, but in a partial-ordered relation induced by the local
density and distance between the data; and the label of a center
can be regarded as the contribution of its followers. Starting from
the assumption, we develop a highly efficient non-iterative label
propagation algorithm based on a novel data structure named
as optimal leading forest (LaPOLeaF). The major weaknesses of
the traditional GSSL are addressed by this study. We further
scale LaPOLeaF to accommodate big data by utilizing block
distance matrix technique, parallel computing, and LocalitySensitive Hashing (LSH). Experiments on large datasets have
shown the promising results of the proposed methods.
Index Terms—Optimal leading forest, semi-supervised learning, label propagation, partial order assumption.

We ponder the possible reasons of these limitations and
argue that the crux is these models treat the relationship among
the neighboring data points as “peer-to-peer”. Because the
data points are considered equal significant to represent their
class, most GSSL objective functions try optimizing on each
data point with equal priority. However, this “peer-to-peer”
relationship is questionable in many situations. For example,
if a data point xc lies at the centering location of the space of
its class, then it will has more representative power than the
other one xd that diverges more from the central location, even
if xc and xd are in the same K-NN or (-NN) neighborhood.
This paper is grounded on the partial-order-relation assumption: the neighboring data points are not in equal status,
and the label of the leader (or parent) is the contribution
of its followers (or children). The assumption is intuitively
reasonable since there is an old saying: “a man is known by
the company he keeps”. The labels of the peripheral data may
change because of the model or parameter selection, but the
labels of the core data are much more stable. Fig. 1 illustrates
this idea.

I. I NTRODUCTION

L

ABELS of data are laborious or expensive to obtain,
while unlabeled data are generated or sampled in tremendous size in big data era. This is the reason why semisupervised learning (SSL) is increasingly drawing the interests
and attention from the machine learning society. Among the
variety of many SSL model streams, Graph-based SSL (GSSL)
has the reputation of being easily understood through visual
representation and is convenient to improve the learning performance by exploiting the corresponding matrix calculation.
Therefore, there has been a lot of research works in this regard,
e.g., [1], [2], [3].
However, the existing GSSL models have two apparent
limitations. One is the models usually need to solve an
optimization problem in an iterative fashion, hence the low
efficiency. The other is that these models have difficulty in
delivering label for a new bunch of data, because the solution
for the unlabeled data is derived specially for the given graph.
With newly included data, the graph has changed and the
whole iterative optimization process is required to run once
again.
Corresponding author: Guoyin Wang.
Ji Xu is with School of Information Science and Technology, Southwest
Jiaotong University. E-mail: alanxuch@hotmail.com
Guoyin Wang is with Chongqing Key Laboratory of Computational Intelligence, Chongqing University of Posts and Telecommunications. E-mail:
wanggy@ieee.org
Manuscript received October 30, 2017.

Unlabeled data
Labeled data

Fig. 1: partial-order-relation assumption: the label of the center
can be regarded as the contribution from the labels of its
followers. Therefore, one can safely infer herein that the left
unlabeled point is a triangle and the right one is a pentagram.
This paper proposes a non-iterative label propagation algorithm taking our previous research work, namely local density
based optimal granulation (LoDOG), as starting point. In
LoDOG, the input data was organized as an optimal number
of subtrees. Every non-center node in the subtrees is led by its
parent to join the microcluster the parent belongs to. In [4],
these subtrees are called Leading Tree. The proposed method,
Label Propagation on Optimal Leading Forest (LaPOLeaF),
performs label propagation on the structure of the relatively
independent subtrees in the forest, rather than on the traditional
nearest neighbor graph.
Therefore, LaPOLeaF exhibits several advantages when
compared with other GSSL methods: (a) the propagation is

IEEE TRANSACTIONS ON XXXX

2

performed on the subtrees, so the edges under consideration
are much more sparse than that of nearest neighbor graph;
(b) the subtrees are relatively independent to each other, so
the massive label propagation computation is easier to be
parallelized when the size of samples is huge;
(c) LaPOLeaF performs label propagation in a non-iterative
fashion, so it is of high efficiency.
Overall, LaPOLeaF algorithm is formulated in a simple way
and the empirical evaluations show the promising accuracy
and very high efficiency.
The rest of the paper is organized as follows. Section II
briefly reviews the related work. The model of LaPOLeaF is
presented in details in Section III. Section IV describes the
method to scale LaPOLeaF for big data. Section V analyzes
the computation complexity and discusses the relationship to
other researches, and Section VI describes the experimental
study. We reach a conclusion in Section VII.

where Z is the regression matrix that describes the relationship
between raw samples and anchors; ΩG (f ) = 21 f T Lf . L =
D − W is the Laplacian matrix.
Wang proposed a hierarchical AGR method to address
the granularity dilemma in AGR , by adding a series of
intermediate granular anchor layer between the finest original
data and the coarsest anchor layer [3].
One can see that the underlying philosophy is still the
two assumptions. Slightly different from the two assumptions,
Ni proposed a novel concept graph harmoniousness, which
integrates the feature learning and label learning into one
framework (framework of Learning by Propagability, FLP)
[2]. The objective function has only one term in FLP, yet
also needs to obtain the local optimal solution by alternately
running iterative optimizing procedure on two variables.

II. R ELATED STUDIES
A. Graph-based semi-supervised learning (GSSL)
Suppose an undirected graph is denoted as G = (V, E, W ),
where V is the set of the vertices, E the set of edges, and
W : E 7→ R is the mapping from an edge to a real number
(usually defined as the similarity between the two ending
points). GSSL takes the input data as the vertices V of the
graph, and places an edge ei,j between two vertices (vi , vj )
if (vi , vj ) are similar or correlated. The basic idea of GSSL
is propagating the labels of the labeled samples Xl to the
unlabeled Xu with the constructed graph. The propagation
strength between vi and vj on each edge is in proportion to
the weight Wi,j .
Almost all the existing GSSL works on two fundamental
assumptions. One is called “clustering assumption”, meaning
that the samples in the same cluster should have the same
labels. Clustering assumption usually is used for the labeled
sample set. The other is called “manifold assumption”, which
means that similar (or neighboring) samples should have similar labels. Manifold assumption is used for both the labeled
and unlabeled data.
Starting from the two assumptions, GSSL usually aims at
optimizing an objective function with two terms. However,
the concrete components in different GSSL models vary. For
example, in [5] the objective function is

X = {x1 , x2 , ..., xN } denotes the dataset . I = {1, 2, ..., N }
is the index set of X. di,j is the distance (under any metric)
between xi and xj .



min
F

l+u
1
1
1X
W ij √ F i − p F j
2 i,j=1
di
dj

2


+µ

l
X

B. optimal leading forest

Definition 1. Local Density. [6] The local density of xi is
di,j 2
P
computed as ρi =
e−( dc ) , where dc is the cut-off
j∈I\{i}

distance or band-width parameter.
Definition 2. leading node and δ-distance. If xli is the nearest
neighbor with higher local density to xi , then xli is called the
leading node of xi . Formally, li = arg min{dij |ρj > ρi },
j

denoted as xli = η(xi ) for short. di,li is called the δ-distance
of xi , or simply δi .
We store all the xli in an array named as LN.
Definition 3. leading tree (LT) [4]. If ρr = max {ρi }, xli =
1≤i≤N

η(xi ) . Let an arrow start from xi , i ∈ I\{r}, and end at xli .
Thus, X and the arrows form a tree T . Each node xi in T
(except xr ) tends to be led by xli to join the same cluster
xli belongs to, unless xi itself makes a center. Such a tree is
called a leading tree.
Definition 4. η operator [7]. For any non-root node x in an
LT, there is a leading node p for x. This mapping is denoted
as η(x) = p.
we denote η( η( ...η(•))) = η n (•) for short.
|
{z
}
n times

||F i − Y i ||

i=1

(1)
where F is the label indication matrix; di is the sum of the
i − th row of W ; Y i is the label of the i − th labeled data.
Liu proposed an Anchor Graph Regulation (AGR) approach
to predict the label for each data point as a locally weighted
average of the labels of anchor points [1]. In AGR, the function
is
c
c
X
1X
||Zl aj − y j ||2 + γ
ΩG (Zaj )
min
Q(A) =
2 j=1
A=[a1 ,...,ac ]
j=1
(2)

Definition 5. partial order in LT [7]. Suppose xi , xj ∈ X ,
we say xi ≺ xj , if f ∃m ∈ N + such that xj = η m (xi ).
Definition 6. center potential. Let γi denote the potential of
xi to be selected as a center, γi is computed as γi = ρi ∗ δi .
Intuitively, if an object xi has a large ρi (means it has
many near neighbors) and a large δi (means relatively far from
another object of larger ρ), then xi would have great chance
to be the center of a collection of data.
Pedrycz proposed the principle of justifiable granularity
indicating that a good information granule (IG) should has
sufficient experiment evidence and specific semantic [8], [9],

IEEE TRANSACTIONS ON XXXX

3

[10]. That is, there should be as many as possible data points
included in an IG, and the closure of the IG should be compact
and tight in geometric perspective.
Following this principle, we have proposed a local density
based optimal granulation (LoDOG) method to build justifiable granules accurately and efficiently [11]. In LoDOG, we
construct the optimal IGs of X by disconnecting the corresponding Leading tree into an optimal number of subtrees. The
optimal number Ngopt is derived via minimizing the objective
function:

points merged in the fat node xi in the subtree, if the node
is derived as an information granule, after some granulation
methods such as local sensitive hashing (LSH) (e.g., [12] [13])
or others. If there is no granulation performed before the LT
construction, all popi are assigned with constant 1.
LaPOLeaF is designed to consist of three stages after the
OLF has been constructed, namely, from children to parent
(C2P), from root to root (R2R), and from parent to children
(P2C). The idea of these stages is illustrated in Figure 2.
x5

L1

min Q(Ng ) = α ∗ H(Ng ) + (1 − α)
α

X

r1

DistCost(Ωi ), (3)

L2

x6

r2

Ng

x3

x4

L1

x7
x8

x9

L2

i=1

where DistCost(Ωi ) =

|ΩP
i |−1

L3

{δj |xj ∈ Ωi \R(Ωi )}.

j=1

Here, Ng is the number of IGs; α is the parameter striking
a balance between the experimental evidence and semantic;
Ωi is the set of points included in ith granule; | • | returns the
cardinality of a set; H(•) is a strictly monotonically increasing
function used to adjust the magnitude of Ng to well match
Ng
P
that of
DistCost(Ωi ). This function can be automatically
i=1

selected from a group of common functions such as logarithm
functions, linear functions, power functions, and exponential
functions; R(Ωi ) is the root of the granule Ωi as a leading
tree.
We used LoDOG to construct the Optimal Leading Forest
(OLF) from the dataset. The readers are referred to [11] for
more details of LoDOG.

x2

x10 x11

x1

(a) From children to parent

(b) From root to root

L3

(c) From parent to children

Fig. 2: Diagram of non-iterative label propagation on the
subtrees of an FNLT. (a) x3 gets its label as the weighted
summation of the L2 and L1 , and the label of x5 is computed
likely in a cascade fashion. (b) r1 is the root of an unlabeled
subtree. In this situation, we have to borrow the label for r1
from r2 . If r2 is not labeled either, this “borrow” operation
will be transitively carried out (see Section III-A2 for details).
(c) After the previous two stages, all roots of the subtrees are
guaranteed being labeled. Then also under the guidance of 4,
all the unlabeled children will get their label information in a
top-down fashion.
To decide the layer number for each node, one can easily
design a hierarchical traverse algorithm (see Appendix) for the
sub-leading-tree.

Definition 7. optimal leading forest (OLF). Ngopt leading
trees can be constructed from the dataset X by using LoDOG
method. All the leading trees are collectively called optimal
leading forest.

A. Three key stages of label propagation in LaPOLeaF

The concept of OLF is used to determine the localized
ranges of label propagation on the whole leading tree of X .
That is, OLF indicates where to stop propagating the label of
a labeled datum to its neighbors.

Definition 8. unlabeled (labeled) node. A node in the subtree
of the OLF is an unlabeled node (or the node is unlabeled),
if its label vector is a zero vector. Otherwise, i.e., if its label
vector has at least one element greater than zero, the node is
called a labeled node (or the node is labeled).

III. L ABEL P ROPAGATION ON O PTIMAL L EADING F OREST
(L A POL EA F)
LaPOLeaF first makes a global optimization to construct the
OLF, then performs label propagation on each of the subtrees.
Following the aforementioned partial relation assumption, the
relationship between the children and their parent is formulated as 4, and each stage of the label propagation of
LaPOLeaF will be guided by this formula.
P
popi
i Wi ∗ Li
, where Wi =
.
(4)
Lp = P
W
dist(i,
p)
i
i
where Lp is the label vector of the parent for a Kclassification problem. Li is the label vector of the i-th child
w.r.t. the current parent. That the k-th element equals one and
all others equal to zero represents a class label of the k-th
class, 1 ≤ k ≤ K. For regression problems, Li and Lp are
simply scalar value. popi is the population of the raw data

1) From children to parent:

Definition 9. unlabeled (labeled) subtree. A subtree in OLF
is called an unlabeled subtree (or the subtree is unlabeled),
if every node in this tree is not labeled. Otherwise, i.e., if the
leading tree contains at least one labeled node, this tree is
called a labeled subtree (or the subtree is labeled).
Since the label of a parent is regarded as the contribution
of its children, the propagation process is required to start
from the bottom of each subtree. The label vector of an
unlabeled children is initialized as vector 0, therefore it does
not contribute to the label of its parent. Once the layer index
of each node is ready, the bottom-up propagation can start to
execute in a parallel fashion for the labeled subtrees.
Proposition 1. After C2P propagation, the root of a labeled
subtree must be labeled.
Proof. According to the definitions of labeled node and
labeled subtree, and the procedure of C2P propagation, a

IEEE TRANSACTIONS ON XXXX

parent is labeled if it has at least one child labeled after the
corresponding round of the propagation. The propagation is
progressing sequentially along the bottom-up direction, and the
root is the parent at the top layer. Therefore, this proposition
obviously holds.
2) From root to root: If the labeled data are rare or unevenly
distributed, there would be some unlabeled subtrees. In such
a case, we must borrow some label information from other
labeled subtrees. Because the label of the root is more stable
than other nodes, the root of an unlabeled subtree ru should
borrow label information from a root of a labeled subtree rl .
However, there must be some requirements for rl . To keep
consistence with our partial order assumption, rl is required
to be superior to ru and is the nearest root to ru . Formally,

4

(a)

(b)

xp
Layeri

x1 x2

x3 x4

Layeri+1

xp'

L p  L p 

W1L1  W2 L2
W1  W2

x3 x4 L3  L4  Lp

Fig. 3: Illustration of the virtual parent idea in P2C propagation
stage. (a) A parent and 4 children (first 2 are labeled, and the
last two are unlabeled). (b) To compute the labels for x3 and
x4 , the labeled nodes (xp , x1 and x2 ) are replaced by a virtual
parent xp0 .

Algorithm 1: LaPOLeaF Algorithm
Input: Dataset X = Xl ∪ Xu
Output: Labels for Xu
1 Part 1: //Preparing the OLF;
rl = arg min{dist(ru , ri )|ru ≺ ri },
(5)
2 Compute distance matrix Dist for X ;
ri ∈RL
3 Compute local density ρ (Definition 1) ;
where RL is the set of labeled roots.
4 Compute leading nodes LN, δ-distance δ (Definition 2) ;
If there exists no such rl for a particular ru , we can conclude
5 Compute representation power γ using (6);
that the root rT of the whole leading tree constructed from X
6 Split the LT into OLF using objective function (3), return
(before splitting into a forest) is not labeled. So, to guarantee
roots RT and nodes set in each subtree;
every unlabeled root can successfully borrow a label, we only
7 Build adjacent List for each subtree;
need to guarantee rT is labeled.
8 Part 2: //Label propagation on the OLF;
If rT is unlabeled after C2P propagation, we consider the
9 Decide the level Number for each node using a
label-borrowing trick for rT as well. However, there is no other
hierarchical traverse approach (See Appendix);
root satisfying rT ≺ rl , so we modify (5) a little to borrow
10 C2P propagation using (4);
T
label from rl for rT :
11 R2R propagation using (5) and (6);
12 P2C propagation using (7);
l
rT = arg min{dist(rT , ri )}
(6) 13 Return Labels for X
u
ri ∈RL
The R2R propagation is executed for the unlabeled roots in
the representation-power-ascending order.
3) From parent to children: After the previous two stages,
all root nodes of the subtrees are labeled. In this P2C propagation, the labels are propagated in a top-down fashion, i.e., the
labels are sequentially propagated from top layer to bottom
layer and this process can be parallelized on the independent
subtrees.
We need consider two situations: a) for a parent xp , all m
children xi , 1 ≤ i ≤ m, are unlabeled. Here, We simply assign
Li =Lp , because this assignment directly satisfies (4) no matter
what value each Wi takes. b) for a parent xp , without loss of
generality, assume the first ml children are labeled, the other
mu children are unlabeled. In this situation, we generate a
virtual parent xp0 to replace the original xp and the ml labeled
children. Using (4), we have
Pml
j=1 Wj ∗ Lj
(7)
Lp0 = Lp − Pml
j=1 Wj

Then, the mu unlabeled children can be assigned with the
label Lp0 like in the first situation. The concept of virtual
parent is illustrated in Fig. 3.
B. LaPOLeaF algorithm
We present the overall algorithm of LaPOLeaF here, including some basic information about OLF construction.

C. An Example on double-moon dataset
We generate the double-moon dataset of 600 data points,
300 in each moon, to illustrate the main stages in LaPOLeaF,
helping the readers to build an intuitive impression of this
method. 5 labeled points are randomly selected in each moon.
In the first step, the OLF with Ngopt = 40 was constructed
using the steps described in Part 1 of Algorithm 1 (Fig.
4a). Here, the parameters in LoDOG are set as: {percent =
0.7, α = 0.5, H(x) = x}. The root of each subtree is marked
by yellow face and red edge. It is easily observed that the
edges appear in the OLF is much sparser than that in other
GSSL methods based on nearest neighbors [2] [5].
In the C2P propagation (Fig. 4b), the nodes in the subtrees
are firstly tagged by layer index. The sub-tree with greatest
height has 14 layers. After the bottom-up label propagation,
the root of each labeled subtree becomes labeled. And other
nodes on the path from the initial labeled node to its corresponding root are labeled as well. There are 44 nodes labeled
now. Yet the unlabeled subtrees remain unchanged.
Fig. 4c shows that in R2R propagation stage, each unlabeled root borrowed the labeled from its nearest neighboring
root with higher density. The green arrows show the label
borrowing information, with arrow head indicating the label
owner.
The P2C propagation can be fully parallelized because all

IEEE TRANSACTIONS ON XXXX

5

Root
Class One
Class Two

(d)

(c)

owner
borrower

D. Deriving the label for a new datum
A salient advantage of LaPOLeaF is that it can obtain the
label for a new datum (let us denote this task as LXNew)
in O(n) time. Because (a) the leading tree structure can
be incrementally updated in O(n) time, and the LoDOG
algorithm can find Ngopt in O(n) time, OLF can be updated
in O(n) time. And (b) the label propagation on the OLF takes
O(n) time.
The interested reader can refer to our previous work [7],
in which we have provided an detailed description of the
algorithm for incrementally updating the fat node leading tree.
And we also provided therein a proof of the correctness of the
algorithm.
IV. S CALABILITY OF L A POL EA F
To scale the model LaPOLeaF for big data context, we
propose two approaches. One uses parallel computing platform
and divide-and-conquer strategy to obtain an exact solution,
and the other is an approximate approach based on LocalitySensitive Hashing (LSH).
A. Divide and conquer approach
The problem confronted has three aspects with the divideand-conquer strategy. (a) Computing the distance matrix with
O(N 2 ) time complexity. (b) The computation of ρi and δi
needs accessing a row of elements in the whole distance
matrix, so ρ and δ for all data also has the complexity

N1

N2

...

NB

D1,1

D1,2

...

D1,B

D2,2

...

D2,B

N2
...

Fig. 4: An illustrative example of LaPOLeaF on the doublemoon dataset. (a)The OLF constructed from the dataset. (b)
C2P propagation. (c) R2R propagation. The green arrows
indicate the borrower and the owner when an unlabeled root
borrows label from another root. All the roots are labeled after
this stage. (d) P2C propagation. The color saturation reflects
the value of the maximal element in a label vector. The closer
to 1 the value is, the higher saturation of the color.

N1

...

Unlabeled Data

...

(b)

(a)

of O(N 2 ). (c) The distances between the centers should be
prepared in advance in the R2R propagation stage, since the
memory of a single computer is not able to accommodate the
whole distance matrix for a large dataset, and the distances
between centers can not be retrieved directly from the whole
distance matrix. Apart from the three parts, other steps in
LaPOLeaF are all linear to N and usually could run on a
single machine.
1) Compute the distance matrix in parallel: Distance matrix puts a considerable burden on both the computation time
and the capacity of memory. As an example, we have computed the required memory capacity for the distance matrix of
100,000 samples is over 37GB, even when a distance is stored
in a float of 4 Bytes.
Here, we propose a divide-and-conquer method for exactly
(not approximately) computing large distance matrix, whose
idea is illustrated in Fig. 5.

...

the roots in OLF are labeled and the propagation within each
subtree is independent from others. Using the discussion in
Section III-A3, all the unlabeled non-root nodes are labeled
(as in Fig. 4d).

NB
(a)

...

DB,B

(b)

Fig. 5: (a) If the whole dataset is divided into two subsets,
then the whole distances between any pair of points can be
computed through tree parts. The first two parts correspond
to the full connections within the two subgraphs respectively
(green and red curves), and the third part corresponds to the
full connections within the complete bipartite graph (black
lines). (b) The number of subsets is generalized from 2 to B.
Note that although we try to balance the size of each subset,
all the Ni are not necessarily equal. therefore, while Di,i is
always square matrix, Di,j for i 6= j may NOT be square
matrix.
Although the computing of distance matrix is of O(N 2 )
complexity, the positive message is that the mainstream CPU
manufacturers (such as Intel and AMD) and scientific computing softwares (such as Matlab and R) have made great
efforts to accelerate the matrix operations. For the distance
metric l2 -norm, instead of computing the distance between
the objects pair by pair, we formulate the distances for the
full connections between the two parts of a bipart graph as
in Theorem 2. For small datasets of 1,000 instances and 8
attributes on Matlab 2014a, matrix computing runs about 20
times faster than pairwise distance computing.
Theorem 2. The Euclidean distance matrix for the full connections within a complete bipartite graph is given by elemente D
e is computed via
wise square root of D.
e = P .2 × 1d×n + 1m×d × Q.2 − 2P × Q> ,
D

(8)

where P and Q are the matrixes formed with data points (as
row vectors) in each part of the bipartite graph, whose sizes

IEEE TRANSACTIONS ON XXXX

6

are m × d and n × d, respectively. P.2 is the element-wise
square for P .
e i,j , we can write
Proof. Considering an element D
e i,j =
D

d
P
k=1

=

d
P

p2i,k +

d
P
k=1

2
qj,k
−2

d
P

pi,k qj,k

k=1

(9)

TABLE I: Complexity comparison.
Methods

Graph construction

Label propagation

LLGC

O(n2 )

O(n3 )

FLP

O(n2 )

O(T1 Kn2 + T2 K 2 n2 )

AGR

O(T mn)

O(m2 n + m3 )

HAGR

O(T mh n)

O(m2h n + m3h )

LaPOLeaF

O(n2 )

O(n)

2

(pi,k − qj,k )

k=1

B. Approximate approach with LSH
Thus, the theorem is proven.
Ideally, the distance matrix of arbitrarily-sized data can be
computed in this way provided there are enough computers.
However, if the computers are not adequate for the dataset at
hand, one can turn to the second approach LSH.
2) Computing ρ, δ and nneigh in parallel: Because of the
additive characteristic of local density ρi , the whole ρ vector
can be computed in a fully parallel fashion on each computing
node, when the whole matrix are split in belt and stored
separately on different computing nodes. Suppose there are
B blocks of the distance matrix for N samples, then we have
ρ=

B
X
b=1

N
X

ρb =

B
X

[ρb1 , · · · , ρbN ],

(10)

b=1

where a is a random vector and b is a random real number
sampled uniformly from the interval (0, r]. For angular similarity, the hash function could be [14]
hv (x) = sgn(v > x)

2

Db (i, j)
) ), i = 1, ..., N/B.
exp(−(
ρbi =
dc
j=1

As mentioned above, if there are no adequate computers
for a given large dataset to run exact LaPOLeaF, then it is
reasonable to merge some closely neighboring data points into
one bucket by employing LSH techniques [12]–[14]. The basic
idea of LSH is that the nearest located neighbors have high
probability to share the same hash code (viz. collide with each
other), and the far away data points are not likely to collide.
For different distance metrics, we need different hash functions. For l2 -norm, the hash function is [12]


a·v+b
(13)
ha,b (v) =
r

(11)

where, ρb is the local density vector of N elements w.r.t.
distance matrix block b, and Db (i, j) is the (i, j) element in
the bth distance matrix block.
Unlike on a single computer, where δ can be computed with the guidance of sorted ρ, computing each
δi in parallel has to access the whole ρ and all
Df loor(i/bSize)+1 (mod(i, bSize), j) for 1 ≤ j ≤ N with
bSize = N/B.
3) Prepare the distance matrix for centers: The R2R propagation stage needs to access the distances between any pair of
centers (roots of the subtrees in OLF), denoted as distCenters .
If the distance matrix is stored on a centralized computer,
distCenters can be extracted directly from the whole distance
matrix D. However, when the D is stored in a distributed
system, it is divided into k blocks, denoted as Db , b = 1, ..., k.
Each Db is stored on a different computing node, and the index
range for instances whose distances are indicated by Db is
[(b−1)∗bSize+1, b∗bSize]. Usually we have bSize = N/k,
except for the last matrix block.
Therefore, to extract distCenters from the distributed D,
one has to sort firstly the centers according to their index in
ascending order, then get the distance entry between Center i
and Center j via
distCenters (i, j) = Df loor(i/bSize)+1 (mod(i, bSize), j);
(12)
By sorting the centers, each distance matrix needs to be
accessed only once to get distCenters .

(14)

where, v is a random vector, and sng(•) is the sign function.
Ji improved the work in [14] by introducing the GramSchmidt orthogonalization process to the random vector group
to form a representing unit named as “super-bit”.
After running LSH algorithm on the original dataset, the
instances will be put into many buckets. Then each bucket is
treated as a fat node by LaPOLeaF, and the number of the
data lying in the ith bucket is the popi in (4).
V. T IME COMPLEXITY AND RELATIONSHIP TO RELATED
WORKS

A. Complexity analysis
By investigating each step in Algorithm 1, we find out that
except the calculation of the distance matrix requires exactly
n(n − 1)/2 basic operations, all other steps in LaPOLeaF has
the linear time complexity to the size of X . When compared
to LLGC [5], FLP [2], AGR [1], and HAGR [3], LaPOLeaF
is much more efficient, as listed in Table I. In Table I, n is the
size of X ; T is the number of iterations; K is the number of
classes; mh is the number of points on the hth layer. Empirical
evaluation in Section VI verified this analysis.
Please note that although we write O(n3 ) for straightforward computation of matrix inverse, the complexity could be
reduced to O(n2.373 ) [15] [16].
It is also worthwhile to compare the efficiency of LXNew.
For LaPOLeaF, LXNew only require a linear time complexity
w.r.t. the size of the existing OLF. However, for the traditional
GSSL methods, LXNew requires the running time as long as
for labeling all data points in X .

IEEE TRANSACTIONS ON XXXX

7

TABLE II: Information of the datasets in the experiments
# Instances

# Attributes

# Classes

Iris

150

4

3

Wine

178

13

3

Yeast

1,484

8

8

MNIST

70,000

784

10

Activity

43,930,257

16

6

B. Relationship discussion
The label propagation in LaPOLeaF is a heuristic algorithm that lacks an optimization objective. Hence it offers no
mathematical guarantee to achieve best solution in this stage.
However, we argue that the optimization has been moved
forward to the OLF construction stage. Since we obtained
an optimal partial ordered structure of the whole dataset, we
believe that an iteration optimization which regards the data
as in peer-to-peer relation is no longer compulsory. This is the
difference between LaPOLeaF and other GSSL methods.
Meanwhile, LaPOLeaF can be regarded as an improved
version of K-NN. In K-NN, the K nearest neighbors is
considered as a spherical-shaped information granule and the
unlabeled data are assigned label with a voting strategy. The
parameter K is set by a user , and the result is quite sensitive to
the choice of K. By contrast, in LaPOLeaF, the information
granules are arbitrarily-shaped leading trees and the size of
each tree is automatically decided by the data and LaPOLeaF,
so usually the sizes are different. Because OLF better captures
the nature of the data distribution, and the label propagation
is reasonably designed, LaPOLeaF constantly outperforms KNN.
VI. E XPERIMENTAL STUDIES
The efficiency and effectiveness of LaPOLeaF is evaluated
on five real world datasets, among which three are small data
from UCI machine learning repository and the other two are
larger scaled. The information of the datasets is shown in
Table II. The 3 small datasets are used to demonstrate the
effectiveness of LaPOLeaF, and the other two datasets are used
to show the scalability of LaPOLeaF through the measures of
parallel computing and Locality-sensitive Hashing (LSH).
The experiment for the small datasets and Activity data is
conducted on a personal computer with an Intel i5-2430M
CPU and 16GB DDR3 memory. MNIST data is learned both
on the PC and a Spark cluster of eight work stations.
A. UCI small datasets
With the 3 small sized UCI datasets, namely, Iris, Wine,
and Yeast, it is shown that LaPOLeaF achieves competitive
accuracy while the efficiency is much higher, when compared
with the classical semi-supervised learning methods Linear
Discriminant Analysis (LDA) [17], Neighborhood Component
Analysis (NCA) [18] , Semi-supervised Discriminant Analysis
(SDA) [19], and Framework of Learning by Propagability
(FLP) [2]. The parameter configuration and some experimental
details, such as the distance metric chosen and preprocessing
method, for all the 5 datasets are listed in Table III.

2.972

x 10

5

α=0.3
α=0.35
α=0.4

2.9715

Objective function value

Dataset

2.971
2.9705

Ng=65
2.97

Ng=153
2.9695
2.969

Ng=405
0

100

200

300

Number of granuels

400

500

Fig. 6: Objective function values v.s. the number of information granules in LoDOG method for MNIST data.
The accuracies for the competing models on the 3 datasets
are shown in Table IV, from which one can see that LaPOLeaF
achieves the best accuracy twice and the accuracy is comparable even when the other method FLP wins on the dataset
Wine.
The main purpose of LaPOLeaF is not to improve the
accuracy of GSSL, but to improve its efficiency by getting
rid of the paradigm of iterative optimization to achieve the
minimum value of an objective function. LaPOLeaF exhibits
very high efficiency, for example, it completes the whole SSL
process within 0.27 seconds for Iris dataset, on the mentioned
personal computer.
B. MNIST dataset
The MNIST dataset contains 70,000 handwriting digits
images in total, about 7,000 for each digit (’0’-’9’). To
emphasize the effectiveness of LaPOLeaF itself, we directly
use the original pixel data as the learning features as in [20].
Since the distance matrix is oversized, we applied the divideand-conquer technology described in Section IV-A. The whole
dataset X is equally divided into 7 subsets, so the size of each
distance matrix block is 10,000×70,000.
After computing the matrix blocks Db , the vector group
{ρ, δ, nneigh}, the OLF of X can be constructed by running
LoDOG algorithm on a single machine. The parameters and
some intermediate results for the two datasets MNIST and
Activity are detailed in the last two rows of Table III. The
objective function values for choosing Ngopt with MNIST data
is shown in Fig. 6.
10 labeled samples are randomly chosen from each digit,
and the accuracies achieved by LaPOLeaF and the stateof-the-art method Hierarchical Anchor Graph Regularization
(HAGR) [20] is listed in Table V.
One can see that LaPOLeaF achieves competitive accuracy
on MNIST data. However, the highlight of LaPOLeaF is its
efficiency. LaPOLeaF complete the whole learning process,
including the OLF construction and the three stages in label
propagation, within 48 minutes on the personal computer. The
time consumptions are detailed in Table VI.
C. Activity dataset
The Activity dataset is from the domain of Human Activity
Recognition [21]. It includes the monitoring data sampled

IEEE TRANSACTIONS ON XXXX

8

TABLE III: Parameters configuration for the 5 datasets
opt

Dataset

percent

α

H(x)

Ng

|Xl |

Preprocessing

Distance

Iris

2

0.25

x + 30

8

6

Z-score

Euclidean

Wine

2

0.4

x1.2

8

6

NCA DR†

Euclidean

Yeast

5

0.1

55 + 2x

7

16

Z-score

Cosine

MNIST

10

0.3

2.97 × 105 + x1.2

405

100

none

Euclidean

Yeast

8

0.3

0.2x

347

16

see Fig. 7

Cosine

†: DR here is the abbreviation of dimensionality reduction.

TABLE IV: Accuracy comparison on the 3 small datasets.
Method

Iris

Wine

Yeast

LDA

66.91±25.29

62.05±

19.00±8.70

NCA

92.28±3.24

83.10±9.70

32.76±6.32

SDA

89.41 ±5.40

90.89±5.39

37.00±6.89

FLP

93.45±3.09

93.13±3.32

40.03±5.40

LaPOLeaF

94.86±4.57

90.68±5.32

42.28±2.36

TABLE V: Accuracies by HAGR and LaPOLeaF on MNIST
Method

HAGRbase

HAGR

LaPOLeaF

Accuracy

79.17±1.39

88.66±1.23

84.92±2.35

by the accelerators and gyroscopes built in the smart phones
and smart watches. Since there are many different models of
phones and watches, and the sampling frequency, accuracy are
therefore different, the data collected are heterogeneous. The
dataset contains 43,930,257 observations and 16 attributes in
each.
1) preprocessing: Because the raw data in separated into
4 comma separated values (.csv) files, we firstly perform
preprocessing as shown in Fig. 7.
Align & merge

ECDF feature

Normalization

SB-LSH

Fig. 7: Preprocessing flow chart of the Activity data.
i) The records from the four .csv files are aligned and
merged into one file. We use the subject ID and equipment
ID to align the data from different files, and the difference in
sampling frequency are dealt with by interpolation.
ii) Since that empirical cumulative distribution function
(ECDF) feature has been reported to outperform FFT and PCA
in HAR task [22]. We compute the ECDF feature from the
original data and use it in the future learning. Conventionally,
the time frame is set to 1 second, and the overlapping ratio is
50%. Since the major sampling frequency is 200, we include
200 observations in one frame to compute ECDF, and move
forward 100 observations when one row of ECDF features
is computed. In this way, the time granularity of activity
recognition is half a second. The segment number is set to
5, so the resulting dimensionality of the feature is 6*5=30.
ECDF has reduced the size of Activity data to 439,302 rows
(about 1% of the original).
TABLE VI: Running time of each stage on MNIST
Stage

Db

ρ

δ

LoDOG

LP$

Time(s)

952

1483

436

15

3

$: LP here is the abbreviation of label propagation.

iii) Use min-max normalization to normalize every column
of the data.
iv) Because of the large number of samples and features,
we employ LSH (specifically, SB-LSH [13]) to tackle the
two problems at the same time. With SB-LSH method, we
empirically set the parameter depth of a super-bit as 5, number
of super-bits as 6. LSH will further reduce the amount of data
rows by merging the collided samples into the same bucket.
For example, the resultant number of hash buckets is 771 for
the subset of subject a and smart phone model Nexus4 1,
compared with the ECDF feature rows’ number of 3,237. The
number of samples contained in each hash bucket is treated
as the weight to compute ρ in OLF construction.
2) LaPOLeaF for Activity data: After the preprocessing
stage, Activity data has been transformed into 84,218 rows of
categorical features. Each row consists of the series of hash
bucket number, and the weight of each row is the number of
ECDF rows that share the same hash bucket number sequence.
The distance matrix is computed based on the ECDF features
in the bucket sequence rather than the hash codes themselves,
and cosine distance is used. The parameters in OLF is set
as (α = 0.3, percent = 8, H(x) = 0.2x). This parameter
configuration leads to the number of subtrees in OLF is 347.
With the constructed OLF, we run C2P, R2R, and P2C
propagation sequentially, taking randomly selected 120 (20 per
class) labeled ECDF features as the labeled data Xl . The final
accuracy achieved by LaPOLeaF is 86.27±3.36.
VII. C ONCLUSIONS
The existing GSSL methods have two weaknesses. One is
low efficiency due to the iterative optimization process, and
the other is inconvenience two predict the label for newly
arrived data. This paper firstly made a sound assumption that
the neighboring data points are not in equal positions, but
lying in a partial-ordered relation; and the label of a center
can be regarded as the contribution of its followers. Based on
this assumption and our previous work named as LoDOG, a
new non-iterative semi-supervised approach called LaPOLeaF
is proposed. LaPOLeaF exhibits two salient advantages: a)
It has much higher efficiency than the sate-of-the-art models
while keep the accuracy comparable. b) It can deliver the labels
for a few newly arrived data in a time complexity of O(N ),
where N is the number of the old data forming the optimal
leading forest (OLF). To enable LaPOLeaF to accommodate
big data, we proposed an exact divide-and-conquer approach
and an approximate locality-sensitive-hashing (LSH) method.
Theoretical analysis and empirical validation have shown the
effectiveness and efficiency of LaPOLeaF. We plan to extend

IEEE TRANSACTIONS ON XXXX

LaPOLeaF in two directions: one is to apply it into the real
world big data mining problem with tight time constraints, and
the other is to improve the accuracy while keeping the high
efficiency unchanged.
A PPENDIX
We provide the algorithm of deciding the layer index of
each node in a tree using a queue data structure.
Algorithm 2: Decide the layer index of nodes in a tree.
Input: The root of T and the adjacent list AL of the tree.
Output: The layer indices for the nodes LayerInd[n].
1 Initialize an empty queue theQue;
2 EnQue(theQue,T);
3 LayerInd[T]=1;
4 while !IsEmpty(theQue) do
5
QueHead=DeQue(theQue);
6
if AL[QueHead]!=NULL then
7
EnQue(AL[QueHead]);
8
LayerInd[AL[QueHead]]= LayerInd[QueHead]+1;
9
end
10 end
11 Return LayerInd[];
The time complexity of Algorithm 2 is O(n), because the
basic operations for each node are EnQue() and DeQue().
ACKNOWLEDGMENT
This work has been supported by the National Key Research and Development Program of China under grants
2016QY01W0200 and 2016YFB1000905, the National Natural Science Foundation of China under grant 61572091.
R EFERENCES
[1] W. Liu, J. He, and S.-F. Chang, “Large graph construction for scalable
semi-supervised learning,” in Proceedings of the 27th international
conference on machine learning (ICML-10), pp. 679–686, 2010.
[2] B. Ni, S. Yan, and A. Kassim, “Learning a propagable graph for semisupervised learning: Classification and regression,” IEEE Transactions on
Knowledge and Data Engineering, vol. 24, no. 1, pp. 114–126, 2012.
[3] M. Wang, W. Fu, S. Hao, H. Liu, and X. Wu, “Learning on big
graph: Label inference and regularization with anchor hierarchy,” IEEE
Transactions on Knowledge and Data Engineering, vol. 29, no. 5,
pp. 1101–1114, 2017.
[4] J. Xu, G. Wang, and W. Deng, “DenPEHC: Density peak based efficient
hierarchical clustering,” Information Sciences, vol. 373, pp. 200–218,
2016.
[5] D. Zhou, O. Bousquet, T. N. Lal, J. Weston, and B. Schölkopf, “Learning
with local and global consistency,” in Advances in neural information
processing systems, pp. 321–328, 2004.
[6] A. Rodriguez and A. Laio, “Clustering by fast search and find of density
peaks,” Science, vol. 344, no. 6191, pp. 1492–1496, 2014.
[7] J. Xu, G. Wang, T. Li, W. Deng, and G. Gou, “Fat node leading tree for
data stream clustering with density peaks,” Knowledge-Based Systems,
vol. 120, pp. 99–117, 2017.
[8] W. Pedrycz and W. Homenda, “Building the fundamentals of granular
computing: A principle of justifiable granularity,” Applied Soft Computing, vol. 13, no. 10, pp. 4209–4218, 2013.
[9] W. Pedrycz, G. Succi, A. Sillitti, and J. Iljazi, “Data description: A
general framework of information granules,” Knowledge-Based Systems,
vol. 80, pp. 98–108, 2015.
[10] X. Zhu, W. Pedrycz, and Z. Li, “Granular data description: Designing
ellipsoidal information granules,” IEEE Transactions on Cybernetics,
DOI: 10.1109/TCYB.2016.2612226, 2016.

9

[11] J. Xu, G. Wang, T. Li, and W. Pedrycz, “Local density-based optimal
granulation and manifold information granule description,” IEEE Transactions on Cybernetics, DOI: 10.1109/TCYB.2017.2750481, 2017.
[12] M. Datar, N. Immorlica, P. Indyk, and V. S. Mirrokni, “Locality-sensitive
hashing scheme based on p-stable distributions,” in Proceedings of the
twentieth annual symposium on Computational geometry, pp. 253–262,
ACM, 2004.
[13] J. Ji, J. Li, S. Yan, B. Zhang, and Q. Tian, “Super-bit localitysensitive hashing,” in Advances in Neural Information Processing Systems, pp. 108–116, 2012.
[14] M. S. Charikar, “Similarity estimation techniques from rounding algorithms,” in Proceedings of the 34th Annual ACM Symposium on Theory
of Computing, pp. 380–388, ACM, May 2002.
[15] V. V. Williams, “Breaking the coppersmith-winograd barrier,” 2011.
[16] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction
to Algorithms, Third Edition. MIT Press, 2009.
[17] P. N. Belhumeur, J. P. Hespanha, and D. J. Kriegman, “Eigenfaces
vs. fisherfaces: recognition using class specific linear projection,” IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 19,
pp. 711–720, Jul 1997.
[18] J. Goldberger, S. Roweis, G. Hinton, and R. Salakhutdinov, “Neighborhood component analysis,” Advances in Neural Information Processing
Systems, pp. 513–520, 2004.
[19] D. Cai, X. He, and J. Han, “Semi-supervised discriminant analysis,” in
2007 IEEE 11th International Conference on Computer Vision, pp. 1–7,
Oct 2007.
[20] M. Wang, W. Fu, S. Hao, H. Liu, and X. Wu, “Learning on big
graph: Label inference and regularization with anchor hierarchy,” IEEE
Transactions on Knowledge and Data Engineering, vol. 29, pp. 1101–
1114, May 2017.
[21] A. Stisen, H. Blunck, S. Bhattacharya, T. S. Prentow, M. B. Kjærgaard,
A. Dey, T. Sonne, and M. M. Jensen, “Smart devices are different:
Assessing and mitigatingmobile sensing heterogeneities for activity
recognition,” in Proceedings of the 13th ACM Conference on Embedded
Networked Sensor Systems, pp. 127–140, ACM, 2015.
[22] N. Y. Hammerla, R. Kirkham, P. Andras, and T. Ploetz, “On preserving
statistical characteristics of accelerometry data using their empirical
cumulative distribution,” in Proceedings of the 2013 International Symposium on Wearable Computers, pp. 65–68, ACM, 2013.

Ji Xu received the B.S. from Beijing Jiaotong
University, Beijing, in 2004, and the M.S. from
Tianjin Normal University, Tianjin, China, in 2008.
Now he is a Ph.D. candidate with Southwest Jiaotong University, Chengdu, China. His research interests include data mining, granular computing and
machine learning. He has published a number of
papers in refereed international journals such as
IEEE Transactions on Cybernetics and Information
Sciences, etc.

Guoyin Wang received the B.S., M.S., and Ph.D.
degrees from Xian Jiaotong University, Xian, China,
in 1992, 1994, and 1996, respectively. He worked at
the University of North Texas, USA, and the University of Regina, Canada, as a visiting scholar during
1998-1999. Since 1996, he has been working at the
Chongqing University of Posts and Telecommunications, where he is currently a professor, the Director
of the Chongqing Key Laboratory of Computational
Intelligence, the Director of the National International Scientific and Technological Cooperation Base
of Big Data Intelligent Computing, and the Dean of the Graduate School. He
was the President of International Rough Sets Society (IRSS) 2014-2017. He
is the Chairman of the Steering Committee of IRSS and the Vice-President of
Chinese Association of Artificial Intelligence. He is the author of 15 books,
the editor of dozens of proceedings of international and national conferences,
and has over 200 reviewed research publications. His research interests include
rough set, granular computing, knowledge technology, data mining, neural
network, and cognitive computing.

