Przemyslaw Spurek
email: przemyslaw.spurek@ii.uj.edu.pl

Jacek Tabor
email: jacek.tabor@ii.uj.edu.pl

Krzysztof Misztal
email: krzysztof.misztal@fis.agh.edu.pl

Przemys Law Spurek

Weighted Approach to Projective Clustering

Keywords: Projective clustering, Karhunen-Loéve Transform, PCA, kmeans

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Projective clustering is a part of the large subspace clustering family, which general aim lies in dividing the given high-dimensional data-set into clusters. For the survey, motivation and further references on subspace clustering we refer to [START_REF] Agarwal | k-Means projective clustering[END_REF][START_REF] Kriegel | Clustering high-dimensional data: A survey on subspace clustering, pattern-based clustering, and correlation clustering[END_REF][START_REF] Parsons | Subspace clustering for high dimensional data: a review[END_REF][START_REF] Vidal | Subspace clustering[END_REF].

A typical application of the projective clustering concerns splitting of a given data-set S into clusters with respect to k affine subspaces. A typical illustrative motivation to create such an algorithm is given by a problem to determine linear components of the data obtained in acoustical experiments [START_REF] Kamisiński | The Study of Sound Scattering Structures for the Purposes of Room Acoustic Enhancement[END_REF] (see Fig. 1(a)). In general, this kind of data contains two linear components: the first connected with sound (product in experiments) which linearly disappears being absorbed by the walls and the air, and the second consists of noise connected with empty space. To use statistical analysis, we have to extract both of them.

The simplest solution, see [START_REF] Grim | Multimodal discrete Karhunen-Loéve expansion[END_REF], lies in a natural modification of k-means: instead of finding k centers which best represent the data, we find k subspaces of given dimension. In other words we change the cost function from the squareddistance from the center of the cluster to the squared distance from the affine space which best represents it and can be found by PCA [START_REF] Jolliffe | Principal component analysis[END_REF]. To explain it graphically let us consider the following example. Fig. 2(a) represents three lines in the plane. The method given in [START_REF] Grim | Multimodal discrete Karhunen-Loéve expansion[END_REF] will split the data into three lines. However, it will not work sufficiently well on data given in 2(b), which we would like to split into a line and a circle.

In this paper we generalize the approach from [START_REF] Grim | Multimodal discrete Karhunen-Loéve expansion[END_REF] to allow the weights over subspaces of different dimensions, which gives a partial solution to the above mentioned problem. The resulted algorithm we call (ω, k)-means, where ω represents the weights, and k the number of clusters we are interested in. In analogy to the case of k-means, we obtain a version of the Voronoi diagram (see next section). In the simplest form our algorithm needs the number of clusters k and the weight parameter ω (for ω = (1, 0, . . . , 0) we obtain the k-means while for k = 1 we obtain the PCA). Detailed investigation of the results of (ω, k)-means allows to determine the optimal dimensions of subspaces, see algorithm in Section 4. The results of this approach we apply to image compression (Section 4, Example 3). In such a case we obtain almost twice the compression of the classical approach.

Generalized Voronoi Diagram

The Voronoi diagram is one of the most useful data structures in computational geometry, with applications in many areas of science [START_REF] Okabe | Spatial tessellations: concepts and applications of Voronoi diagrams[END_REF]. For the convenience of the reader and to establish the notation we shortly describe the classical version of the Voronoi diagram (for more details see [START_REF] Klein | Concrete and Abstract Voronoi Diagrams[END_REF]). For N ∈ N consider R N with the standard Euclidean distance and let S be a finite subset of R N . For p, q ∈ S such, that p = q, let

B(p, q) = {x ∈ R N : p -x = q -x }, (1)
D(p, q) = {x ∈ R N : p -x < q -x }. (2)
The hyperplane B(p, q) divides R N into two subsets, one containing points which are closer to point p then q (D(p, q)), and the second one containing points which are closer to point q then p (D(q, p)) -see Figure 3 Now we proceed to the description of our modification of the Voronoi diagram. We divide the space R N with respect to affine subspaces of R N . For n ≤ N let

E n (R N ) := {(v 0 , . . . , v n ) ∈ (R N ) n+1 : v i , v j are orthonormal for i, j > 0, i = j}.
Thus v 0 denotes a center of affine space we consider, while v 1 , . . . , v n is the orthonormal base of its "vector part". From the geometrical point of view the

element v = (v 0 , v 1 , . . . , v n ) ∈ E n (R N ) represents the affine space v 0 + lin(v 1 , . . . , v n ) = aff(v 0 , v 1 , . . . , v n ).
We modify equations ( 1) and ( 2), by using distance between a point and affine subspace.

Definition 1. Let n < N and let v ∈ E n (R N ), ω = (ω 0 , . . . , ω n ) ∈ [0, 1] n+1 such that n j=0 ω j = 1 be given. For x ∈ R N let DIST ω (x; v) :=   n j=0 ω j dist(x; aff(v 0 , . . . , v j )) 2   1/2 , (3)
where dist(x; V ) denotes the distance of the point x from the space V .

In formula [START_REF] Fisher | On grouping for maximum homogeneity[END_REF], ω = (ω 0 , . . . , ω n ) is interpreted as vector of weights, where ω k denotes the weight of the affine subspace of dimension k.

Remark 1. It is easy to notice, that DIST has the following properties:

for v ∈ E n (R N ) and ω = (0, . . . , 0, 1) ∈ [0, 1] n+1 we obtain that DIST ω (x; v) is a distance between the point x and affine space aff(v); if v 0 = 0 and ω = (0, . . . , 0, 1) then DIST ω is a distance between point and linear space generated by (v 1 , . . . , v n ); if ω = (1, 0, . . . , 0) then DIST ω is the classical distance between x and v 0 :

DIST ω (x; v) = x -v 0 . -if ω =     0, . . . , 0 k , 1 l -k , . . . , 1 l -k k-l , 0, . . . , 0     for k < l, then DIST ω describes
the mean distance between x and subspaces of dimension from k to l.

Corollary 1. Formula (3) can be reformulated as follows To optimize calculations we define

(DIST ω (x; v)) 2 = n j=0 ω j x -v 0 2 - j i=1 x -v 0 ; v i 2 = n j=0 ω j x -v 0 2 - n j=0 ω j j i=1 x -v 0 ; v i 2 . (a) (b) (c) (d) (e)
v1 = x -v 0 ; v 1 2 , vj = vj-1 + x -v 0 ; v j 2 ,
and since ω j = 1 we can simplify our formula to

(DIST ω (x; v)) 2 = x -v 0 2 - n j=0 ω j vj .
Now we are ready to define our generalization of the Voronoi diagram. Let S be a finite subset of E n (R N ) and ω ∈ [0, 1] n+1 , ω j = 1, where n ≤ N . For p, q ∈ S such, that p = q, let B ω (p, q) := {z ∈ R N : DIST ω (z; p) = DIST ω (z; q)}, D ω (p, q) := {z ∈ R N : DIST ω (z; p) < DIST ω (z; q)}.

The set B ω (p, q) divides the space R N into two subsets, containing points which are closer to p then to q (D ω (p, q)) and points which are closer to q then p (D ω (q, p)).

Definition 2. Let n ∈ N, n < N be fixed. Let S be a finite subset of E n (R N ) and ω ∈ [0, 1] n+1 , n j=0 ω j = 1 be given. For p ∈ S the set D ω (p, S) := q∈S : q =p D ω (p, q)
of all points that are closer to p than to any other element of S is called the (open) generalized Voronoi region of p with respect to S.

Applying this definition we obtain a new type of Voronoi diagram. Figure 4 presents a generalized diagram on the plane for different weights changing from ω = (1, 0) to ω = (0, 1). In general we obtain that the boundary sets usually are not polygons but zeros of quadratic polynomials. The same happens in R 3 even for ω = (0, 1) see the Fig. 5, where we show points with equal distance from two lines.

Generalization of the k-means method

Clustering is a classical problem of the division of a set S ⊂ R N into separate clusters, or in other words, into sets showing given type of behavior.

One of the most popular and basic method of clustering is the k-means algorithm. By this approach we want to divide S into k clusters S 1 , . . . , S k with minimal energy. For convenience of the reader and to establish the notation we shortly present the k-means method.

For a cluster S and r ∈ R N we define the cost function

E(S, r) := s∈S s -r 2
which we interprets as an energy. We say that the point r best "describes" the set S if energy is minimal, more precisely, if

E(S, r) = inf r∈R N {E(S, r)}.
It is easy to show that barycenter (mean) of S minimizes the function E(S, •) (for more information see [START_REF] Ding | K-means clustering via principal component analysis[END_REF][START_REF] Fisher | On grouping for maximum homogeneity[END_REF]). The above consideration can be precisely formulated as follows:

Theorem 1 (k-means). Let S be a finite subset of R N . We have

E(S, µ(S)) = inf r∈R N {E(S, r)}
where µ(S) := 1 cardS s∈S s denotes the barycentre of S. Thus in the k-means the goal is to find such clustering S = S 1 ∪ . . . ∪ S k that the function E(S 1 , . . . , S k ) = k j=1 E(S j , µ(S j )) is minimal. In this paper we consider generalization of k-means algorithm similar to that from the previous section concerning the Voronoi diagram. Instead of looking for points which best "describe" clusters we seek n dimensional subspaces of R N .

Let

S ⊂ R N and ω ∈ [0, 1] n+1 , ω j = 1 be fixed. For v ∈ E n (R N ) let E ω (S, v) := s∈S DIST 2 ω (s, v).
We interpret the function E ω (S, v) as an energy of the set S respectively to the subspace generated by v. If the energy is zero, the set S is subset of affine space generated by v. We say that v best "describes" the set S if the energy is minimal, more precisely if

E ω (S, v) = inf v∈En(R N )
{E ω (S, v)}.

To obtain an optimal base we use a classical Karhunen-Loeve transform (called also Principal Component Analysis, shortly PCA), see [START_REF] Jolliffe | Principal component analysis[END_REF]. The basic idea behind the PCA is to find the coordinate system in which the first few coordinates give us a "largest" possible information about our data.

Theorem 2 (PCA). Let S = {s 1 , . . . , s m } be a finite subset of R N . Let

M(S) := (v 0 , . . . , v N ) ∈ E N (R N )
be such that

-v 0 = µ(S); -v 1 , . . . , v N are pairwise orthogonal eigenvectors of [s 1 -v 0 , . . . , s m -v 0 ] • [s 1 - v 0 , . . . , s m -v 0 ]
T arranged in descending order with respect to the eigenvalues.

For every n < N and ω ∈ [0, 1] n+1 we have

E ω (S, M k (S)) = inf v∈En(R N ) {E ω (S, v)},
where M k (S) := (v 0 , . . . , v k ).

Thus given ω ∈ [0, 1] n+1 , ω j = 1, in (w, k)-means our goal is to find such clustering S = S 1 ∪ . . . ∪ S k that the function

E ω (S 1 , . . . , S k ) := k j=1 E ω (S j , M n (S)) (4)
is minimal. Consequently (ω, k)-means algorithm can be described as follows:

stop condition choose ε > 0 initial conditions choose randomly points {s 1 , . . . , s k } ⊂ S obtain first clustering (S 1 , . . . , S k ) by matching each of the points s ∈ S to the cluster such that s -s j 2 is minimal repeat let E = E ω (S 1 , . . . , S k ) compute vectors v 1 , . . . , v k , which best "describe" clusters, by the PCA method (v j = M n (S j ))

obtain new clustering (S 1 , . . . , S k ) by adding each of the point s ∈ S to the cluster such that DIST ω (s, v j ) is minimal until E -E ω (S 1 , . . . , S k ) < ε As is the case in the classical k-means, our algorithm guarantees a decrease in each iteration but does not guarantee that the result will be optimal. It is easy to notice that the above method has following properties:

1. for ω = (1, 0, . . . , 0) we obtain the classical k-means, 2. for n = 1 we get Karhunen-Loéve transform.

Example 1. As already mentioned in Section 2, the k-means do not find a global minimum and strongly depends on initial selection of clusters. In our case, this effect can be even more visible. Consider the case of circle C in R 2 with 4 clusters and ω = (0, 1). The picture, see Figure 6

Determining the dimensions of subspaces

In this section we present the method of determining the optimal dimensions subspaces which describe clusters. By using various value of parameter ω we are able to compress data by remembering different coordinates in each cluster.

Let k (the number of clusters) and ω (weight parameter) be fixed. As a result of the (ω, k)-means algorithm for the dataset S we obtain k clusters {S 1 , . . . , S k } and k coordinate systems {v 1 , . . . ,

v k } ⊂ E n (R N ). For v ∈ E n (R N ) and n 0 ≤ n we define sub-base of dimension n 0 by v n0 = (v 0 , . . . , v n0 ).
We choose n 1 , . . . , n k ∈ {1, . . . , N } and we compress the data of S by replacing each element s ∈ S i by its orthogonal projection on a suitable subspace spanned on v ni .

Let S 1 , . . . , S k and {v 1 , . . . , v k } ⊂ E n (R N ) be a result of the (ω, k)-mean algorithm. For parameters n 1 , . . . , n k we consider the compression error Comp err(n 1 , . . . , n k ) :=

k i=1 s∈Si dist 2 (s; v i ni ) 1/2
.

Let ε > 0 be given a maximal allowed error. We want to find the minimal number of parameters to "compress" the data with compression error below ε.

Observe that if we approximate S i by its projection onto subspaces of dimension n i , the total number of parameters is given by

n 1 • card(S 1 ) + . . . + n k • card(S k ).
The procedure of determining respective dimensions n 1 , . . . , n k of clusters, such that Comp err(n 1 , . . . , n k ) < ε can be formulated as follows:

1. Apply the (ω, k)-means algorithm with given k (in general this parameter should be chosen respectively to data structure) and ω (which describe possible dimensions of clusters). 2. In each cluster S 1 , . . . , S k determinate eigenvalue of covariance matrix λ j 1 ≥ . . . ≥ λ j N for j = 1, . . . , k.

Put

Λ := λ l1 1 , . . . , λ l kN kN .

4. Sort the eigenvalues increasingly

Λ (•) = λ l (1)
(1) , . . . , λ l (kn) (kn) . 5. Let n := sup n :

n i=1 λ l (i) (i) • m l (i) ≤ ε ,
where m i = card(S i ), for i = {1, . . . , k}. 6. We define n 1 , . . . , n k by

n j = card λ l ( i) (i) : such that l (i) = j and (i) > n .
Before we show that this algorithm gives good accuracy we present following theorem Lemma 1 ( [START_REF] Jolliffe | Principal component analysis[END_REF]). Let S = {x 1 , . . . , x n } be subset of R N . By {λ 1 , . . . , λ n } we denote eigenvalues corresponding to eigenvectors {v 1 , . . . , v n } of matrix cov([x 1 , . . . , x n ]).

Then

x∈S dist 2 (x, v k ) = n i=k+1 λ i • n, where v k = {v 0 , v 1 , . . . , v k }.

Now by simple calculations we have

Comp err(n 1 , . . . , n k ) = k i=1 s∈Si dist(s; v ni ) 2 1/2 =   k i mi j=ni+1 λ i j • m i   1/2 = n i=1 λ l (i) (i) • m l (i) 1/2 < ε.
Example 2. Consider the dataset containing points group around a segment (200 points) and circle (200 points) -see Figure 7. In first step we fix ε = 2.87 (which gives 5.5% of total error3 ). Then to start our algorithm we have to choose the parameters k and ω. In our example we want to obtain two clusters, so we fix k = 2. Moreover the first cluster should represent the one-dimension data and the second two-dimension. So we put ω = (0, 1 2 , 1 2 ). Outcome obtained at the end of calculation is presented in Table 2. Cluster S 1 corresponds to points grouped along interval, and the S 2 -along circle. Now by steps 3-6 we have

Λ (•) = {0.001, 0.001, 0.040, 1.885, 2.124, 9.241} , n = 3, n 1 = 1, n 2 = 2.
Consequently, we get 2 • 199 parameters for S 2 and 1 • 201 for S 1 .

At the end of this section we presen our algorithm on the example of the classical Lena picture. Since in picture compression we expect the data to have lower dimensional structure 4 , we narrow our consideration to subspaces of dimension between 10-20 by choosing, according to Remark 1, By applying points 3-6 we obtain:

-1 • 2375 parameters for the first cluster, -2 • 151 parameters for the second cluster, -5 • 880 parameters for the third cluster, -4 • 229 parameters for the fourth cluster, -3 • 461 parameters for the fifth cluster.

As we see by use our method we have to remember 9376 parameters. If we fix n 1 = . . . = n 5 such that Comp err(n 1 , . . . , n 5 ) < ε (4 first eigenvalues for each cluster) we obtain 16384 parameters -which is all most twice as much as in our method.

Sample implementation of (ω, k)-means algorithm prepared in Java programming language is available at [START_REF] Misztal | Implementation of the (ω, k)-means algorithm[END_REF].

Fig. 1 .

 1 Fig. 1. Linear component in acoustical data structure. Fig. 1(a) -original data. Fig. 1(b) -outcome from (ω, k)-means algorithm for k = 2, ω = (0, 1). We extract two linear components in data (black dots match clusters centers with the corresponding lines describing those clusters, vertical line separate sound and background noiseafter 4.3 s).

Fig. 2 .

 2 Fig. 2. Our goal is to create algorithm which will split Fig. 2(b) as two groups of oneand two-dimensional points.

Fig. 3 .

 3 Fig. 3. Graphical presentation of B(p, q), D(p, q) and D(p, S) in R 2 .

  (a). The set D(p, S) := q∈S : q =p D(p, q) of all points that are closer to p than to any other element of S is called the (open)Voronoi region [7] of p with respect to S. Set D(p, S) for N = 2 is the interior of a convex, possibly unbounded polygon (Figure 3(b)). The points on the contour of D(p, S) are those that have more than one nearest neighbour in S, one of which is p. The union V (S) := ∂D(p, S) of all region boundaries is called the Voronoi diagram [7] of S. The common boundary of two Voronoi regions is a Voronoi edge. Two edges meet at a Voronoi vertex.

Fig. 4 .

 4 Fig. 4. Generalized Voronoi diagram for clustering of 3 clusters for different weight vectors Fig. 4(a), ω = (1, 0); Fig. 4(b), ω = ( 3 4 , 1 4 ); Fig. 4(c), ω = ( 1 2 , 1 2 ); Fig. 4(d), ω = ( 1 4 , 3 4 ); Fig. 4(e), ω = (0, 1).

Fig. 5 .

 5 Fig. 5. Generalized Voronoi diagram for ω = (0, 1) and two lines.

Fig. 6 .

 6 Fig. 6. Circle clustering in R 2 for 4 clusters with ω = (0, 1). (ω, k)-means method strongly dependents on initial conditions.

  (a), shows clustering obtained by use (ω, k)-means algorithm. Of course it is a local minimum of E ω , however as we see at Figure 6(b) it is far from being the global minimum.

Fig. 7 .

 7 Fig. 7. Clustering with (ω, k)-means for ω = (0, 1 2 , 1 2 ).

Table 1 .

 1 Outcome of the (ω, k)-means in the case of data from Example 2.Example 3. Let us consider the classical Lena image. First, we interpret photo as a matrix. We do this by dividing a picture into disjoint squares (8 by 8) pixels, where each of them is described (in RGB) by using 3 parameters. Consequently we have vectors from R 192 . Let ε = 427 (which gives 1% of total error) be fixed. We use k = 5. Then we have to choose ω. If we do not have any intuition about possible dimension of cluster we can put

				S1				S2
	µ(Si)	(-0.048, -0.027, 0.0) (-0.004, 0.002, -0.012)
	eigenvector	 	0.0 0.0 -0.845 -0.534 0.534 -0.845	 	 	0.692 -0.722 0.001 0.722 0.692 0.002	 
			1.0	0.0	0.0		0.003	0.0	-1.0
	eigenvalue (9.241, 0.040, 0.001) (2.124,1.885,0.001)
			ω =	1 192	, . . . ,	1 192	.

By total error we understand error obtained in the worst case of compression, when we replace each element in each cluster by barycenter of all data.

That why the compression based on Karhunen-Loéve transform or JPG format gives good results.


Robert Burduk
email: robert.burduk@pwr.wroc.pl

Pawel Trajdos

Construction of Sequential Classifier Using Confusion Matrix

Keywords: Multistage classifier, sequential classifier, confusion matrix

This paper presents the problem of building the decision scheme in the multistage pattern recognition task. This task can be presented using a decision tree. This decision tree is built in the learning phase of classification. This paper proposes a split criterion based on the analysis of the confusion matrix. Specifically, we propose the division associated with an incorrect classification. The obtained results were verified on the data sets form UCI Machine Learning Repository and one real-life data set of the computer-aided medical diagnosis.

Introduction

The classification task may be divided according to its complexity. There are two groups here: one-step and multistage approach. In the one-step approach there is no division into smaller tasks classification. However, a multistage (sequential) approach breaks up a complex decision into a collection of several simpler decisions [START_REF] Ko Lakowska | Fisher Sequential Classifiers[END_REF][START_REF] Mui | Automated classification of nucleated blood cells using a binary tree classifier[END_REF][START_REF] Podolak | Hierarchical classifier with overlapping class groups[END_REF][START_REF] Safavian | A survey of decision tree classifier methodology[END_REF]. Many algorithm build a tree structure in the learning process [START_REF] Penar | Experiments on classifiers obtained via decision tree induction methods with different attribute acquisition cost limit[END_REF][START_REF] Quinlan | Induction on Decision Tree[END_REF]. In other approaches the decision tree structure is fixed before the learning process [START_REF] Manwani | Geometric decision tree[END_REF]. Generally, the synthesis of the multistage classifier is a complex problem. It involves a specification of the following components [START_REF] Safavian | A survey of decision tree classifier methodology[END_REF][START_REF] Kurzyński | Decision Rules for a Hierarchical Classifier[END_REF]:

design of a decision tree structure [START_REF] Woźniak | A hybrid decision tree training method using data streams[END_REF], selection of features used at each non-terminal node of the decision tree [START_REF] Choraś | Image feature extraction methods for ear biometrics-A survey[END_REF][START_REF] Choraś | Content-based retrieval using color, texture, and shape information[END_REF][START_REF] Guyon | An introduction to variable and feature selection[END_REF][START_REF] Rejer | Genetic Algorithms in EEG Feature Selection for the Classification of Movements of the Left and Right Hand[END_REF], the choice of decision rules for performing the classification [START_REF] Burduk | Decomposition of classification task with selection of classifiers on the medical diagnosis example[END_REF].

In particular, this paper discuses a way to design a decision tree structure. The split criterion is based on the confusion matrix. The potential division of the node is associated with the analysis of misclassification in the learning process. In the experiment decision rules are chosen arbitrarily in the entire tree.

The content of the work is as follows. Section 2 introduces the idea of the hierarchical (sequential) classifier. In Section 3 we describe the proposed split criterion. In the next section we present the results of the experiments verified on data sets form UCI repository and one real-life data set of the computer-aided medical diagnosis. The last section concludes the paper.

Hierarchical classifier

The hierarchical classifier contains a sequence of actions [START_REF] Burduk | Classification error in Bayes multistage recognition task with fuzzy observations[END_REF][START_REF] Kurzyński | On the Multistage Bayes Classifier[END_REF]. These actions are simple classification tasks executed in the individual nodes of the decision tree. Some specific features are measured on every non-leaf node of the decision tree. At the first nonleaf node features x 0 are measured, at the second features x 1 are considered and so on. Every set of features comes from the whole vector of features. In every node of the decision tree the classification is executed according to the specific rule. The decisions i 0 , i 1 , ..., i N are the results of recognition in the suitable node of the tree. The design of a decision tree structure is based on split criterion.

In our task of classification the number of classes is equal to N C. The terminal nodes are labeled with the number of the classes from M = 1, 2, ..., N C, where M is the set of labels classes. The non-terminal nodes are labeled by numbers of 0, NC+1, NC+2 reserving 0 for the root-node. Let us introduce the notation for the received model of the multistage recognition [START_REF] Kurzyński | Decision Rules for a Hierarchical Classifier[END_REF]:

-M -the set of internal (nonleaf) nodes, -M i -the set of class labels attainable from the i-th node (i ∈ M), -M i -the set of nodes of the immediate descendant node i (i ∈ M), m i -the node of the direct predecessor of the i-th node (i = 0).

In each interior node the recognition algorithm is used. It maps observation subspace to the set of the immediate descendant nodes of the i-th node [START_REF] Berger | Statistical Decision Theory and Bayesian Analysis[END_REF][START_REF] Duda | Pattern Classification[END_REF]:

Ψ i : X i → M i , i ∈ M . (1)
This approach minimizes the misclassification rate for the particular nodes of a tree. The decision rules at each node are mutually independent. In experiment the decision rules are chosen arbitrarily in the entire tree. Each of the classifiers used in the nodes of the tree takes a decision based on the full set of attributes available in the training set.

In our method of induction, the classification tree is a regular binary tree. This means that on each of the tree nodes there is a leaf or a node has two children.

Induction of the decision tree is performed by the top-down method. This means that it is initiated by the classifier located in the root of the tree. Using the proposed criterion the decision is made whether to continue the division. The process is repeated for the subsequent child nodes of the tree, until the state wherein the nodes in the tree can no longer be divided.

Split criteria

The division of the internal node will be made on the basis of the multidimensional confusion matrix. Specifically, we propose the division associated with an incorrect classification. This division is binary, which means that the node that will meet the criterion of the split will have two child nodes. One of them represents a new internal decision tree node. The second one represents the label of a class that met the appropriate condition. The internal node of the decision tree is analysed in detail via the multidimensional confusion matrix. The columns of the confusion matrix correspond to the predicted labels (decisions made by the classifier in the internal node). The rows correspond to the true class labels. In this matrix the diagonal elements represent the overall performance of each label. The off-diagonal elements represent the errors related to each label. Now we present the split criterion. For every class labels from internal node we create the L × L dimensional confusion matrix. Now we calculate the set of factors W (k l ), where l = 1, 2, ..., L is the number of class labels, according to the formula:

W (k l ) = L m=1,m =l w l,m + L m=1,m =l w m,l . (2)
The example of the confusion matrix is presented in Tab. 1. The division of the node occurs when max l∈1,2,...,L

W (k l ) -min l∈1,2,...,L W (k l ) > T, (3)
where T ∈ [0, 1] is a fixed threshold value. At the threshold one does not make a division of node. If we set this value at the beginning of the experiment, it indicates that the classification process is performed in the one-stage approach.

Experiments

In the experiential research 10 data sets were tested. Nine data sets come from UCI Machine Learning Repository [19]. The tenth comes from the Surgical Clinic Wroclaw Medical Academy and describes the acute abdominal pain diagnosis problem. A set of all the available features was used for all data sets, however, for the acute abdominal pain data set the selection of features has been made in accordance with the suggestions from another work on the topic [START_REF] Burduk | Different decision tree induction strategies for a medical decision problem[END_REF][START_REF] Kurzyński | Diagnosis of acute abdominal pain using three-stage classifier[END_REF]. The numbers of attributes, classes and available examples of the investigated data sets are presented in Tab. 2. Tab. 3 presents the mean error for 5 -N N (5-nearest neighbor) classifier for the selected values of the parameter T . In Tab. 4 we presented the average ranks for all experiments. The average ranks are calculated on the basis of the Friedman test [START_REF] Trawiński | Nonparametric statistical analysis for multiple comparison of machine learning regression algorithms[END_REF]. Each column in this table is attributed to one test. This means that the one-step approach T = 1 was compared in order to propose in the work the sequential approach. In Tab. 4 the lowest average ranks in each group are shown in bold.

Table 3. Avarage error for 5-NN classifier Data set T=0,02 T=0,04 T=0,02 T=0,4 T=0,5 T=1 Acute 0,161 0,161 0,163 0,165 0,168 0,163 Breast 0,428 0,419 0,438 0,489 0,457 0,412 Ecoli 0,128 0,131 0,124 0,131 0,131 0,131 Glass 0,336 0,344 0,333 0,327 0,34 0,343 Irys 0,03 0,035 0,03 0,03 0,035 0,035 Lung 0,472 0,525 0,537 0,458 0,487 0,557 Seeds 0,112 0,112 0,112 0,112 0,112 0,112 Vert. 0,184 0,181 0,178 0,178 0,181 0,184 Wine 0,302 0,303 0,297 0,322 0,325 0,319 Yeast 0,431 0,437 0,44 0,457 0,444 0,443

All classifiers from the group NN (Nearest Neighbor), with properly chosen values T , may improve the quality of classification. The value of this improvement, however, is not significant from the statistical point of view. For the posthoc Bonferroni-Dunn [START_REF] Demšar | Statistical comparisons of classifiers over multiple data sets[END_REF] test the critical difference (CD) for the 27 values of the parameter T and 10 data sets is equal CD = 11, 7. This CD is calculated at α = 0.05. Although the differences of the average rank do not exceed CD, the obtained results can be considered promising because they are close to this value (in particular 5-NN).

Conclusions

In the paper we propose a split criterion based on analizing the confusion matrix. Specifically, we propose the division associated with an incorrect classification. This criterion is used in the design of a decision tree structure in the multistage classifier. With a fulfilled criteria a binary split of the analyzed decision node is carried out. If we set T = 1, it indicates that the classification process is performed in the one-stages approach. Which means that it does not creat a decision tree.

Experiments done in the work show that we have obtain promising results. The proposed approach improves the quality of classification for k-N N group of classifiers. In some cases, the difference of mean ranks obtained by Friedman test is close to the critical difference. In the future work we can use different division criteria in order to design a sequential classifier. For example, the separable linearization [START_REF] Bobrowski | Separable Linearization of Learning Sets by Ranked Layer of Radial Binary Classifiers[END_REF] or MacArthur's [START_REF] Macarthur | On the relative abundance of bird species[END_REF] overlapping niches model can be used in the split criterium.

Table 1 .

 1 The confusion matrix for the nonleaf node i

			estimated
		k1 k2 . . . kL
	k1 w1,1 w1,2 . . . w1,L
	true k2 w2,1 w2,2 . . . w2,L
	. . .	. . .	. . .	. . .

kL wL,1 wL,2 . . . wL,L

Table 2 .

 2 Description of data sets selected for the experiments

	Data set	example attribute class
	Acute Abdominal Pain 476	31	8
	Breast Tissue	106	10	6
	Ecoli	336	7	8
	Glass Identification	214	10	6
	Irys	150	4	3
	Lung Cancer	31	52	3
	Seeds	210	7	3
	Vertebral Column	310	6	3
	Wine	178	13	3
	Yeast	1484	8	10

Table 4 .

 4 Average ranks from Friedman test

	T	3-NN	5-NN	7-NN 9-NN	SVM
	0	15,45	14,1	19	15,85	15,85
	0.02 13.95	7.8	15.1	16.9	15.85
	0.04 9.95	13.15	17	15.8	15.85
	0.06 14.2	13.05	16.7	16.4	15.85
	0.08 11.85	11.15	14	11.6	15.85
	0.1 12.1	13.45	15.15 13.95	14.5
	0.12 16.85	10.95	19.75 11.2	14.5
	0.14 12.65	19.1	11.7	7.85	14.9
	0.16 12.25	12	14.35 16.15	13.2
	0.18 15.75	14.25	15.5	15.1	13.2
	0.2 15.4	8.55	11.7	10.95	13.2
	0.22 12.85	14.7	7.45	16.15	13.2
	0.24 11.05	13.55	16.15 13.65	13.2
	0.26 18.1	16.5	16.1	14.85	14.55
	0.28 16	12.7	14.3	9.85	14.55
	0.3 13.55	13.8	15.35 16.5	14.55
	0.32 12	12.85	11.5	14.15	13.7
	0.34 14	12.75	13.7	16.4	13.15
	0.36 15.2	15.35	13.2	15.9	13.15
	0.38 10.7	13.55	14.3	14.5	13.15
	0.4 16.05	15.75	14.95 9.65	13.15
	0.42 13.15	13.55	14.7	14.35	13.15
	0.44 12.9	18	10.15 12.75	13.15
	0.46 16.6	15.5	12.75 14.1	13.15
	0.48 13.3	17.95	11.1	13.65	13.15
	0.5 16.15	16.95	11.7	14.45	13.15
	1	16	17	10.65 15.35	13.15

Acknowledgments. The work was supported in part by the statutory funds of the Department of Systems and Computer Networks, Wroclaw University of Technology and by the by The Polish National Science Centre under the grant N N519 650440 which is being realized in years 2011-2014.


Lukáš Vojáček
email: lukas.vojacek@vsb.cz

Jiří Dvorský
email: jiri.dvorsky@vsb.cz

Growing Neural Gas -A Parallel Approach

Keywords: growing neural gas, high-dimensional dataset, high performance computing

The paper deals with the high dimensional data clustering problem. One possible way to cluster this kind of data is based on Artificial Neural Networks (ANN) such as SOM or Growing Neural Gas (GNG). The learning phase of the ANN, which is time-consuming especially for large high-dimensional datasets, is the main drawback of this approach to data clustering. The parallel modification, Growing Neural Gas, and its implementation on the HPC cluster is presented in the paper. Some experimental results are also presented.

Introduction

Recently, the issue of high-dimensional data clustering has arisen together with the development of information and communication technologies which support growing opportunities to process large data collections. High dimensional data collections are commonly available in areas like medicine, biology, information retrieval, web analysis, social network analysis, image processing, financial transaction analysis and many others.

Two main challenges should be solved to process high-dimensional data collections. One of the problems is the fast growth of computational complexity with respect to growing data dimensionality. The second one is specific similarity mea-surement in a high-dimensional space. As presented in [START_REF] Beyer | When is "nearest neighbor" meaningful? In: Database Theory[END_REF], Beyer et al. for any point in a high-dimensional space the expected distance, computed by Euclidean measure to the closest and to the farthest point, shrinks with growing dimensionality. These two reasons reduce the effectiveness of clustering algorithms on the above-mentioned high-dimensional data collections in many real applications.

The authors propose an effective data clustering algorithm which is based on Growing Neural Gas (GNG) [START_REF] Martinetz | A "neural-gas" network learns topologies[END_REF]. The computational complexity is resolved by the parallel implementation of GNG. Some technical problems have to be resolved in order to effectively train such kind of neural network using an High Performance Computing (HPC) cluster with MPI. The traditional serial approach to training GNG is also considered in the paper. The serial learning GNG algorithm is used for benchmarking the parallel version of GNG. In other words, parallel GNG has to produce the same network and should be an order of magnitude faster in the ideal case.

Artificial Neural Networks

Self Organizing Maps

Self Organizing Maps (SOM), also known as Kohonen maps, were proposed by Teuvo Kohonen in 1982 [START_REF] Kohonen | Self-Organization and Associative Memory[END_REF]. SOM is a kind of artificial neural network that is trained by unsupervised learning. Using SOM, the input space of training samples can be represented in a lower-dimensional (often two-dimensional) space [START_REF] Kohonen | Self Organizing Maps[END_REF], called map. SSuch model is efficient in structure visualization due to its feature of topological preservation using a neighbourhood function.

SOM consists of two layers of neurons (see Fig. 1): an input layer that receives and transmits the input information and an output layer, the map that represents the output characteristics. The output layer is commonly organized as a two-dimensional rectangular grid of nodes, where each node corresponds to one neuron. Both layers are feed-forward connected. Each neuron in the input layer is connected to each neuron in the output layer. A real number, or weight, is assigned to each of these connections.

Growing Neural Gas

The principle of this neural network is an undirected graph which need not be continuous. Generally, there are no restrictions on the topology. The graph is generated and continuously updated by competitive Hebbian Learning [START_REF] Martinetz | Competitive hebbian learning rule forms perfectly topology preserving maps[END_REF][START_REF] Prudent | An incremental growing neural gas learns topologies[END_REF]. According to the pre-set conditions, new neurons are automatically added and connections between neurons are subject to time and can be removed. GNG can be used for vector quantization by finding the code-vectors in clusters [START_REF] Holmström | Growing Neural Gas Experiments with GNG, GNG with Utility and Supervised GNG[END_REF], biologically influenced [START_REF] Sledge | Growing neural gas for temporal clustering[END_REF], image compression, disease diagnosis. GNG works by modifying the graph, where the operations are the addition and removal of neurons and edges between neurons. An example of the operation is shown in Figure 2 To understand the functioning of GNG, it is necessary to define the algorithm. The algorithm described by Algorithm 1 is based on the original algorithm [START_REF] Fritzke | A growing neural gas network learns topologies[END_REF] [3], but it is modified for better continuity in the SOM algorithm. The description of the algorithm has been divided for convenience into two parts. In the first part of Algorithm 1 the overall functionality is described. The second part of Algorithm 2 describes one iteration, which means the descriptions of variables without dependence on time.

Remark The notation used in the paper is briefly listed in Table 1.

Algorithm 1 Growing Neural Gas algorithm 3. Perform the one learning iteration according to the algorithm described in Algorithm 2. 4. Reduce error value ei for all neurons Ni using factor β. 5. Returns to step 2, until all input data vector have been used. 6. If t < T return to step 2.

Parallelization

Parallelization focuses on the equally distribution of neurons, where new neurons are allocated to the process with the lowest number of neurons. The advantage of this distribution is constant workload processes. The disadvantage is increase communication between processes.

After analysing the GNG learning algorithm we identified the one most processor time-consuming area. This part was selected as a candidate for the possible parallelization. The selected area is:

Finding BMU -this part of GNG learning can be significantly accelerated by dividing the GNG output layer into smaller pieces. Each piece is then assigned to an individual computation process. The calculation of the Euclidean distance among the individual input vector and all the weight vectors

ec 1 = ec 1 + wc 1 -x 2 (1)
3. Update the weight vector wc 1 of neuron Nc 1

wc 1 = wc 1 + lc 1 (x -wc 1 ) (2)
4. For all neurons N k where exists edge e c 1 k (Nc 1 neighbourhood) (a) Update the weights w k using lnc 1 learning factor

w k = w k + lnc 1 (x -wk) (3) (b) Increase age a kc 1 of edge e c 1 k a kc 1 = a kc 1 + 1 (4) 5.
If there is no edge between neurons Nc 1 and Nc 2 , then create such edge. If the edge exists, the age is set to 0. 6. If any edge has reached the age of amax, it is removed. 7. If there is a neuron without connection to any edge, the neuron is then removed. to find BMU in a given part of the GNG output layer is the crucial point of this part of GNG learning. Each process finds its own, partial, BMU in its part of the GNG output layer. Each partial BMU is then compared with other BMUs obtained by other processes. Information about the BMU of the whole network is then transmitted to all the processes to perform the updates of the BMU neighbourhood.

A detailed description of our approach to the parallelization process is described in Fig. 3.

The parallelization of GNG learning was performed on an HPC cluster, using Message Passing Interface (MPI) technology. MPI technology is based on effective communication between processes. That means that one application can run on many cores. The application uses MPI processes which run on individual cores. The processes are able to send and receive messages and data, communicate etc. Detailed information about HPC and MPI technology is provided, for example, in [8]3 .

Description of Proposed Approach

This subsection includes a detailed description of our approach. Initially, each process reads the training input vectors and on the first process are create two neurons. During the calculation, new neurons are equally distributed on the processes. We do not store the complete GNG graph, we only store parts of it in the appropriate computational nodes.

This approach results in the following significant advantage: neurons are equally distributed among the processes corresponding to the used cores, con- In the main computational phase, one BMU and a second BMU are founded for each input vector. Thus, each processor needs to compute its local BMU and second BMU within its neurons, after which, each local BMU and second BMU (and their position in the network) are shifted onto one process using the MPI function GatherFlattened to determine the global BMU and second BMU. It is possible to use another MPI functions as well, which can provide this selection at one time, but after testing we have found that the experiments took much more time than our presented approach. A global winning BMU and second BMU are then distributed using the MPI function Broadcast on all the processes. Now if there is a edge between the first BMU and the second BMU then age is set to zero otherwise creates edge between this two neurons. Next, the neighbourhood of the BMU in each process is known and, consequently, the weights of the neurons matching are actualized. If the condition is met for adding a new neuron, the process with the lowest number of neurons add a new neuron. A detailed description of adding a new neuron can be found in Fig. 4. This procedure is repeated until all the input training vectors are exhausted (until we have finished one epoch). Two datasets were used in the experiments. The first dataset was commonly used in Information Retrieval -Medlars. The second one was the test data for the elementary benchmark for clustering algorithms [START_REF] Ultsch | Clustering with SOM: U*C[END_REF].

Medlars Dataset

The Medlars dataset consisted of 1,033 English abstracts from a medical science4 . The 8,567 distinct terms were extracted from the Medlars dataset. Each term represents a potential dimension in the input vector space. The term's level of significance (weight) in a particular document represents a value of the component of the input vector. Finally, the input vector space has a dimension of 8,707, and 1,033 input vectors were extracted from the dataset.

Clustering dataset Three training data collections called TwoDiamonds, Lsun and Hepta from the Fundamental Clustering Problems Suite (FCPS) are used.

A short description of the selected dataset used in our experiments is given in Table 2.

First Part of The Experiment

The first part of the experiment was oriented towards a comparison of the standard GNG algorithm and parallel approach to this GNG learning algorithm. The Medlars dataset was used for the experiment. A parallel version of the learning algorithm was run using 2, 8, 16, 24 and 32 MPI processes. The records with an asterisk (*) represents the results for only one process i.e. this is the original serial learning algorithm and there is no network communication. GNG parameters are the same for all experiments and are as follows γ = 200, e w = 0.05, e n = 0.006, α = 0.5, β = 0.0005, a max = 88, M = 1000, δ = 100. The achieved computing time is presented in Table 3. As we can see from Table 3, the computing time depends on the number of used cores as well. With a growing number of processors, the computation effectiveness increases, and the computational time is sufficiently reduced.

Second Part of The Experiment

The second part of the experiments was oriented towards comparing the results obtained by the parallel and standard GNG algorithm. The Clustering dataset was used for the experiment. The parallel version of the learning algorithm was run using 16 MPI processes.

GNG parameters are similar to the previous experiment. There are two changes M = 500 and δ = 25.

Conclusion

In this paper the parallel implementation of the GNG neural network algorithm is presented. The achieved speed-up was very good and the results from the standard and parallel version of GNG are same. So we can say that the operation of the parallel version is correct. However, the effectiveness of a parallel solution is dependent on the division of the output layer. An improper division may cause the communication between processes to be very time consuming.

In future work we intend to focus on the sparse date, use combinations of neural networks for improved result and improved acceleration.

Fig. 1 .

 1 Fig. 1. Basic Schema of SOM

8 .Fig. 2 .

 82 Fig. 2. GNG examples

Fig. 3 .

 3 Fig. 3. Parallel Algorithm

Fig. 4 .

 4 Fig. 4. Parallel Algorithm for adding a neuron

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Results of dataset Target

Table 1 .

 1 Notation used in the paper One iteration of the Growing Neural Gas algorithm 1. Find neurons BMUs neurons Nc 1 and Nc 2 . 2. Update the local error ec 1 of neuron Nc 1

	Symbol Description
	M	Number of input vectors
	n	Dimension of input vectors, number of input neurons, dimension of weight
		vectors in GNG output layer neurons
	N	Current number of neurons in GNG output layer
	Nc 1	The first Best Matching Unit (BM U1), winner of learning competition
	Nc 2	The second Best Matching Unit (BM U2), the second best matching neu-
		ron in learning competition
	wc 1 (t) Weight vector of BM U1
	wc 1 (t) Weight vector of BM U2
	lc 1	Learning factor of BM U1
	lnc 1	Learning factor of BM U1 neighbours
	ei	Local error of output neuron Ni, i = 1, 2, . . . , N
	α	Error ei reduction factor
	β	Neuron error reduction factor
	γ	Interval of input patterns to add a new neuron
	amax Maximum edges age
	aij	Age of edge eij

Nmax Maximum allowed number of neurons in GNG output layer ni i-th input neuron, i = 1, 2, . . . , n Ni i-th output neuron, i = 1, 2, . . . , N eij edge between neurons Ni and Nj for some i, j = 1, . . . , N , where i = j. E set of all edges in GNG G undirected graph describing topology of GNG, G({N1, . . . , NN }, E) T Number of epochs t Current epoch, t = 1, 2, . . . , T X Set of input vectors, X ⊂ R n x(t) Current input vector in epoch t, arbitrarily selected vector from set X x(t) ∈ X, x(t) = (x1, x2, . . . , xn) wk(t) Weight vector of neuron N k , k = 1, 2, . . . , N wk(t) ∈ R n , wk(t) = (w 1k , w 2k , . . . , w nk )

Table 2 .

 2 Fundamental Clustering Problems Suite -selected datasets Experimental Hardware All the experiments were performed on a Windows HPC server 2008 with 6 computing nodes, where each node had 8 processors with 12 GB of memory. The processors in nodes were Intel Xeon 2.27GHz. The topology with the connection between the head node and computing nodes can be found on the web 5 (topology number four). The Enterprise and Private Networks link speed was 1Gbps, the Application link speed was 20Gbps.

	Name	Cases #Vars #Clusters Main Clustering Problem
	Target	770	2	6	outlying clusters
	Lsun	400	2	3	different variances in clusters
	TwoDiamonds	800	2	2	touching clusters

Table 3 .

 3 Computing Time with Respect to Number of Cores, Standard GNG Algorithm, Dataset Medlars

	Cores Computing Time [hh:mm:ss]
	1*	00:35:41
	2	00:17:50
	8	00:04:47
	12	00:03:56
	16	00:03:09
	24	00:02:45
	32	00:02:32

A specification of MPI is available on the web: http://www.mpi-forum.org/ The functions used in the experiment with a link to MPI.Net are available on the web: http://osl.iu.edu/research/mpi.net

The collection can be downloaded from ftp://ftp.cs.cornell.edu/pub/smart. The total size of the dataset is approximately 1.03 MB.

http://technet.microsoft.com/en-us/library/cc719008(v=ws.10).aspx

Acknowledgement This work was supported by the European Regional Development Fund in the IT4Innovations Centre of Excellence project (CZ.1.05/1.1.00/02.0070) and by the Development of human resources in research and development of latest soft computing methods and their application in practice project, reg. no. CZ.1.07/2.3.00/20.0072 funded by Operational Programme Education for Competitiveness, co-financed by ESF and state budget of the Czech Republic.


Anna Bartkowiak

Radoslaw Zimroz

Efficacy of some primary discriminant functions in diagnosing planetary gearboxes

Keywords: condition monitoring, discriminant analysis, healthy state, faulty state, kernels, nonlinearity, gearbox diagnostics

We consider the efficiency of some primary discriminant functions applied in planetary gearbox diagnostics. Real data for planetary gearboxes mounted in bucket wheel excavators working in field condition are elaborated. The aim is to perform condition monitoring (faulty or healthy) of such devices. The raw recorded data (vibration series emitted by the device) were first segmented and transformed to frequency domain using power spectra densities (PSD). Next, 15 variables denoting amplitudes of derived spectra were extracted. This yielded two data matrices A and B of size 1232 × 15, and 951 × 15, representing the faulty and the healthy device appropriately. The data are non-Gaussian and the covariances in both groups differ significantly. Now, using Fisher's discriminant criterion and the kernel methodology, we construct from a learning sample (counting only 600 items) a discriminant function able to provide a discriminant score for distinguishing between the healthy and faulty state of a gearbox. The function proved to be very effective: Both for the learning and the testing data samples (600 and 1483 data vectors respectively) we got 100% correct assignments to the 'faulty' and 'healthy' class, with a conspicuous margin between the two classes. The results are visualized in a 2D plane.

Introduction

Condition monitoring of the state of a gearbox may be done in a multitude of ways, see the survey papers by [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF] or [START_REF] Zimroz | Two simple multivariate procedures for monitoring planetary gearboxes in non-stationary operating conditions[END_REF] for further references. Gearboxes appear worldwide in many devices and machines. The state of a gearbox is important both for economy and safety reasons. A particular attention is paid to planetary gearboxes which represent a complicated architecture of compounds and are costly in exploitation.

The monitoring of the state of such devices is usually carried out by an analysis of vibration recorded by set of sensors (accelerometers) and specialized data acquisition systems (Bruel and Kjaer Pulse system was used here). Additionally, auxiliary channel equipped in tachometer (optical probe used often as rotating shaft speed indicator) was used to perform signal segmentation.

Generally the analysis is performed in two phases: Firstly the vibration signals are preprocessed and set of diagnostic features is extracted. The aim of this phase is to obtain from the recorded vibration data a kind of (multivariate) numerical data data suitable for further analysis. Features extraction in frequency domain (Fourier analysis like the PSD method) is the basic approach for rotating machinery at this stage. Another method often uses for this purpose wavelet transform, however this is done rather for localized damage than for the considered here distributed form of change of condition. The second phase uses various pattern recognition methods, like descriptive statistical methods, neural networks, discriminant functions and classification assignments. Reduction of dimensionality of the data is an important issue in both phases.

A very important point in the analysis is to have proper data for the analysis. For heavy-duty and high power planetary gearboxes the data are difficult to acquire. Therefore quite a lot of researchers use for their analysis data acquired in laboratory conditions using simulations. The behavior of devices working in field conditions might be different from that observed in a laboratory.

Our analysis is based on real data recorded in field conditions. The main contribution is that using these data and applying Kernel Discriminant function (KDF) it is possible to designate a decision boundary allowing for a 100% separation between the faulty and healthy data. The scores of the KDF are visualized in a 2-D plane yielding a visual proof that the two sets of data are really separable.

In the following, Section 2 introduces 3 primary discriminant functions, among them the Fisher's KDF (Kernel Discriminant Function). Section 3 contains a short description of some essential features of the data serving as the basis for our analysis, also some details on constructing the learning and testing samples. In Section 4 results of our analysis, when applying the 3 discriminant function, are presented. Section 5 contains some discussion and closing remarks.

Some primary discriminant functions

There is a multitude of discriminant functions based on diversified criterions. We consider only the two group case when the data are subdivided into two classes (like: 'faulty' or 'healthy') represented by two groups of the observed data. The sought discriminant function should produce for each data vector a score (called also decision score) allowing to assign that vector to one of the considered classes. We will consider only some basic discriminant functions which -to our opinion -are advised to start with, when attempting to build a discriminant function providing a decision boundary between two groups of observed data. In next subsections we consider in turn:

(i) Linear discriminant function formulated and solved as ordinary multivariate regression function with predicted values y = +1 for one of the groups and y = -1 for the other group. (ii) Classical canonical discriminant function based on Fisher's criterion (called often FDA, Fisher's discriminant analysis). (iii) Discriminant functions using kernels providing a nonlinear approach.

Other methods of discriminant analysis -and there are many of them -may be found, among others, in [START_REF] Duda | Pattern Classification[END_REF][START_REF] Hastie | The Elements of Statistical Learning; Data Mining, Inference and Prediction[END_REF][START_REF] Trendafilov | Simple and interpretable discrimination[END_REF].

Denotations:

⋄ X = {x ij }, i = 1, . . . , n, j = 1, 2, . . . , d -the recorded matrix of observations, with d columns denoting variables, and n, the number of rows of the matrix. Each row of X constitutes a data vector x i containing data recorded for the ith individual, called also data item. The vector x i is written as the column vector x i = [x i1 , x i2 , . . . , x id ] T ; it may be also considered as point located in the data space R d . Each of the data vectors {x i } may belong to only one of the two groups of data.

⋄ y = {y i }, i = 1, . . . , n -the vector of group labels. Each data vector has its label indicating to which group (class) the given (ith) data vector is belonging. The group membership is crisp: each data vector may belong only to one of the two considered classes.

Linear discriminant function via linear regression analysis

We consider here the regression model

y i = b 1 x i1 + b 2 x i2 + . . . + b d x id + b 0 + e i , i = 1, . . . , n (1)
The vector y = [y 1 , y 2 , . . . , y n ] T denotes here the predicted variable, which takes specifically only two values: -1, if the ith data vector comes from the faulty class, and +1, if the ith data vector comes from the healthy class.

The coefficients b 1 , . . . , b d , b 0 appearing in (1) are called regression coefficients, they are usually unknown and have to be estimated in a supervised way, that is from a learning sample, which contains known values od the data matrix X and of the vector y. The regression coefficients are commonly estimated using the principle of Least Squares (LS ), which minimizes the following criterion:

SS e = n ∑ i=1 (e i ) 2 = n ∑ i=1 (y i -b 1 x i1 -. . . -b d x id -b 0 ) 2 .
To obtain the solution, that is, the vector

b = [b 1 , b 2 , . . . , b d , b 0 ] T which
yields the minimum of the quadratic form SS e , we have firstly to define X + , the augmented data matrix X, as:

X + = [X, 1 n ],
where 1 n denotes a column of n ones.

With known X + and y, the estimates of b, denoted in the following as b, are obtained as the solution of the following set of linear equations:

(X + ) b = y, (2)
In Matlab the solution ŷ may be obtained simply via the matlab left divide function.

To obtain ŷ, the predicted scores of class assignment of vectors contained in a data matrix say X learn we add to that matrix a column of ones 1 n obtaining the augmented data matrix X + learn , and perform the multiplication ŷlearn = (X + learn ) b.

Similarly, using the vector b yielded by the learning sample, we may obtain the predicted values ŷtest for the test sample

X + test as ŷtest = (X + test ) b.
Predicted values ŷ near to -1 (generally: negative values) indicate similarity to the class 'faulty', and values of ŷ near to +1 (generally: positive values) indicate similarity to the class 'healthy'.

FDA, classical canonical discriminant analysis based on Fisher's criterion

The BW-projections using the between and within groups scatter.

The main idea of the method is to find a vector u = (u 1 , ..., u d ) T serving for construction of a new variable (feature) Z with values z i obtained from the observed data vector x i in the following manner:

z i = x T i u. i = 1, . . . , n. ( 3
)
The values of the new variable, denoted as z i , i = 1, . . . , n, may be considered as projections of the data vectors x i onto a new axis called Z.

In the following we will consider the case when the the entire set of data vectors is subdivided into G = 2 groups (classes). For the ith data point of the gth group (i = 1, . . . , n g , g = 1, 2), we denote the projection of the data point x (g) i as z (g) i , and the general transformation may be read as:

z (g) i = (x (g) i ) T u.
Thus, we obtain G=2 groups of projections:

z (1) 1 , ..., z (1) n1 1st group , z (2) 1 , ..., z (2) n2 2nd group (4)
For values of the variable Z appearing in formula (4) above we may calculate group means z. (g) , (g = 1, 2), and the overall mean z. (.) ≡ z (this is the mean of all data vectors, without subdivision into the 2 groups). Moreover, basing on the values given in (4), we may define and calculate ZW , the within group scatter of the variable Z, and ZB, the between group scatter of the variable Z, as

ZW = 2 ∑ g=1 ng ∑ i=1 (z (g) i -z. (g) ) 2 , ZB = 2 ∑ g=1 n g (z. (g) -z) 2 .
Which projections are the best for the purpose of discriminant analysis? Consider the following issues. For a good discriminant function ⋄ points-projections belonging to different groups of data should be separated as much as possible, ⋄ points-projections belonging to the same groups should be possibly much concentrated around their group-means.

Taking this into account, R.A. Fisher has formulated the following criterion, which is called today Fisher's criterion:

J(z) = ZB ZW = ∑ 2 g=1 n g (z. (g) -z. (.) ) 2 ∑ 2 g=1 ∑ ng i=1 (z (g) i -z. (g) ) 2 . ( 5
)
To indicate for some differentiation between the groups, the criterion J(z), evaluated for real data, should have a value larger than 1; the larger the better. To find the best vector u maximizing the criterion J(z) we express (5) in terms of u and X; usually X is centered to zero, that is, all column means are equal zero. After some algebra manipulations we obtain

J(z) = J(u; X) = u T Bu u T Wu , ( 6
)
where

B d×d = 2 ∑ g=1 n g (x. (g) -x. (.) )(x. (g) -x. (.) ) T , ( 7
)
W d×d = 2 ∑ g=1 ng ∑ i=1 (x l i (g) -x. (g) )(x (g) i -x. (g) ) T . ( 8
)
The above equations show clearly that the criterion J may be considered as a function of the sought unknown vector u and the known data matrix X, which serves for calculations of the between and within group scatter matrices B and W, each of size d × d. Applying Lagrange's methodology and using the additional restraint u T Wu = 1, one arrives to the two-matrices eigenvalue problem:

(B -λW)u = 0. ( 9
)
which is well known in matrix algebra. It is shown, that for G = 2 groups there is an unique vector u satisfying eq. ( 9). Moreover, the constant λ is equal to the quotient J(z) appearing in formula ( 5), evaluated at the solution u. In other words, the constant λ is equal to the ratio of the between group scatter divided by the within group scatter defined above. It satisfies the inequality:

λ ≥ 0
Large values of λ (≫ 1.0) indicate a big differentiation between projections of the two groups of data.

Using the obtained vector u for the projection z = x T u we obtain a new variable Z satisfying Fisher's criterion [START_REF] Bartkowiak | Outliers analysis and one class classification approach for planetary gearbox diagnosis[END_REF]. The technique of constructing such variables is called Canonical Discriminant Analysis (CDA), the obtained projections Z are called Canonical variates or CDA scores. The CDA scores may be calculated and displayed both for the learning and testing samples of data.

The BT projections using the between groups and total scatter. The criterion K(z) formulated defined in formula ( 5) might be formulated in another way. Let T denote the total data scatter matrix:

T d×d = 2 ∑ g=1 ng ∑ i=1 (x l i (g) -x. (.) )(x (g) i -x. (.) ) T , ( 10
)
with x. (.) denoting the overall mean calculated from all the data. Now our aim is to find the vector v maximizing the following index denoting the ratio of the between group to the total scatter

η = v ′ Bv v ′ Tv (11)
The solution of the problem ( 11) is obtained in a similar way as that of (5); it leads to the two-matrix problem

(B -ηT)v = 0. (12)
Now the value eta evaluated at the solution v of ( 12) takes values from the interval 0 ≤ η ≤ 1.0 Big values of η indicate for a good separability of the groups.

Similarly as for Fisher's projections using the B and W matrices, the values z i of the BT-projections are obtained as:

z i = x T i v, ( 13
)
where v is the solution of [START_REF] Liu | Fault level diagnosis for planetary gearboxes using hybrid kernel feature selection and kernel Fisher discriminant analysis[END_REF]. The projections z i obtained that way may be visualized similarly as those obtained from the BW-projections.

Let us emphasize: The canonical discriminant variables represent projections to lower dimension subspaces. They serve in first place for visualization of the projected points. The differentiation of the projected points follows from the assumed Fisher's criterion. The group assignment is relatively easy only in the 2-group case. Having G > 2 groups, one should use for class assignments a specialistic tool, e.g. Mahalanobis distances or k-nearest neighbor methods.

Discriminant functions using kernels

Kernel methods may be defined and used in a number of ways, see for example [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF][START_REF] Mika | Constructing descriptive and discriminative nonlinear features: Rayleigh coefficients in kernel features spaces[END_REF][START_REF] Műller | An introduction to kernel-based learning algorithms[END_REF][START_REF] Baudat | Generalized discriminant analysis using a kernel approach[END_REF][START_REF] Camastra | Kernel Methods for Computer Vision[END_REF][START_REF] Camastra | Kernel Methods for unsupervided Learning[END_REF]. One might say shortly that this is canonical discriminant analysis carried out in an extended space F obtained by a non-linear mapping of the original data. The applied mapping takes into account various nonlinear relations between the observed variables which makes that in the extended feature space F the classical algorithms (like CDA) become more powerful.

Let φ denote the transformation carrying out the mapping of the observed vector x (x ∈ R d ) to the extended feature space F. The function φ is called the mapping function. The extended space F must have defined an operator called scalar product for every pair of its elements φ(x i ), and φ(x j ). The scalar product is defined via so called dot product for pairs of data points obtained as images of the mapping φ of points x i and x j from R d to φ(x i ), and φ(x j ) in F:

k(x i , x j ) = (φ(x i ) • φ(x j )), ( 14
)
where the symbol '•' means the scalar product between the vectors-images φ(x i ) and φ(x j ) evaluated in F. Let us point out, that the kernel matrix K is of size n × n, with n being the size of the learning sample, thus its evaluation may be computationally demanding.

Performing the mapping φ(x) might be quite cumbersome -this depends from the applied mapping φ, in particular from the dimension of F which is usually much larger than that of R d . However, by looking at the transformations carried out using particular mapping functions φ it was found that the elements of the kernel matrix K defined in ( 14) may be evaluated directly from the points x i and x j located in R d without considering their images in F. Then, if we are able to formulate the solution of our task in terms using only the elements of the kernel matrix K, we are the winners: we do not need anymore the mapped points from F. This possibility is called 'the kernel trick'.

For example, this is true when considering as φ the Gaussian or polynomial kernels and having as task canonical discriminant analysis or principal component analysis. Detailed algorithms how to do it in the case of discriminant analysis for G groups of data are shown, e.g., in [START_REF] Műller | An introduction to kernel-based learning algorithms[END_REF] (when using the scatter matrices B, W) or [START_REF] Baudat | Generalized discriminant analysis using a kernel approach[END_REF] (when using the scatter matrices B, T). The algorithms are somehow lengthy and we do not show them here.

More details on the methods and examples of applications may be found in [START_REF] Baudat | Generalized discriminant analysis using a kernel approach[END_REF][START_REF] Camastra | Kernel Methods for Computer Vision[END_REF][START_REF] Camastra | Kernel Methods for unsupervided Learning[END_REF][START_REF] Bartkowiak | Visualizing Some Multi-Class Erosion Data Using Kernel Methods[END_REF][START_REF] Bartkowiak | Visualizing of some multi-class erosion data using GDA and supervised SOM[END_REF][START_REF] Bartkowiak | Visualization of Five Erosion Risk Classes using Kernel Discriminants[END_REF].

Experiments with the gearbox data

In the following we will show an analysis conducted using true data, that is recorded from machines working in field conditions. The data were recorded by Bartelmus and Zimroz [START_REF] Bartelmus | A new feature for monitoring the condition of gearboxes in nonstationary operating systems[END_REF]. We had 2 groups of data contained in 2 data matrices A of size 1232 × 15 and B of size 951 × 15 obtained from vibration data of two have been assigned to the healthy group. This makes the error rate equal to 0.0050. In the testing sample we have 11 incorrect assignments for a total of 1583 items: 10 faulty items were assigned as healthy, and 1 healthy item was assigned as faulty. This makes the error rate equal to 0.0069. Thus the linear discriminant function is a good choice for using it as a discriminant function, especially as it is very easy to carry out.

Results from canonical discriminant analysis

The main result from CDA -when using BW projections -is shown in Figure 2. The second CV2 has discriminative power equal to 0, however is helpful in perceiving the spread of the data points in the two groups.

We have calculated canonical discriminant functions CV1 and CV2 using the BW-projections based on Fisher's Criterion (5) and using own software programmed in Matlab. The main result is shown in Figure 2. Looking at that figure one can grasp at once how different is the spread of points in the two groups of data.

The first canonical variate CV1 has discriminative power lambda=6.8485. We have constructed a second canonical variate which has discriminative power equal to 0, however helps in visualization of the two data clouds representing the faulty and healthy data points.

Concerning class assignments, we have the following results: Learning sample: 3 erroneous assignments in the faulty class. Thus the error rate equals 3/600 = 0.005. Testing sample: There are 24 erroneous assignments in the faulty class. Thus the error rate equals 24/1583 = 0.0154.

Results of Kernel canonical discriminant analysis

The calculations were performed using BT projections elaborated by Baudat and Anouar [START_REF] Baudat | Generalized discriminant analysis using a kernel approach[END_REF] and using their software gda downloaded from www.kernel-machines.org. As previously, the entire observed data were subdivided into the learning sample and test sample -the same, as used in subsections 4.1 and 4.2. The original d=15 dimensional data were mapped to the feature space F where the kernel matrix K = {k(x i , x j }, i, j = 1, . . . , n was constructed from the data vectors belonging to R d using the Gaussian kernels with parameter σ = 0.7:

k(x i , x j ) = exp( -(x i -x j ) T (x i -x j ))/σ (15)
In such a way a non-linear mapping of the original data into F was obtained.

The proper analysis was carried out in the extended space F. The vector v performing the BT-projections onto new axis Z positioned now in F was obtained using the kernel matrix K whose elements were calculated using formula [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. Next, using only information from the kernel matrix K, the projection vector v was constructed. This allowed us to construct new variable Z satisfying the criterion [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF], however valid for the mapped data in F. Next, using the functions spreadGda and plotGda (see [START_REF] Baudat | Generalized discriminant analysis using a kernel approach[END_REF]), it was possible to recover the projections z i , i = 1, . . . , n for each individual data vector from the learning sample (n learn=600 ) and test sample (n test = 1283). The obtained scores z i obtained now are analogous to the scores z i defined for the classical canonical discriminant analysis (see eq. ( 13)) however the obtained now (from spreadgda and plotgda) are projections from the non-linear mapping valid in F.

Discussion and Concluding Remarks

Using real data gathered in field conditions we have investigated the usefulness of some primary discriminant functions in diagnosis of the faulty or healthy state of two planetary gearboxes. Thanks to applying a nonlinear transformation of the gathered data (using Gaussian radial basis functions) to an extended feature space F we have obtained a perfect decision boundary between the faulty and healthy data. The function designating decision boundary was obtained from a relative small learning sample counting 300 faulty and 300 healthy data vectors. The behavior of the obtained decision boundary was checked using a larger test sample counting 1583 data vectors. All of them were correctly classified to the proper group of their origin. We were also able to represent the projections from the extended feature space in an ordinary 2-D plane, as shown in figs. 1 and3.

The elaborated data were not-normal and with covariance matrices differing substantially [START_REF] Zimroz | Investigation on spectral structure of gearbox vibration signals by principal component analysis for condition monitoring purposes[END_REF][START_REF] Zimroz | Multidimensional data analysis for condition monitoring: features selection and data classification[END_REF]. Despite that the constructed discriminants have worked surprisingly effectively.

Fig. 1 .

 1 Fig. 1. Discriminant scores obtained via linear regression method with target values y = -1 for the faulty and y = +1 for the healthy state of the machine. Left panel: learning sample, first 300 faulty and next 300 healthy items. Right panel: testing sample, first 932 faulty and next 651 healthy items.

Fig. 2 .

 2 Fig. 2. Display of two canonical discriminant functions constructed from d=15 dimensional data. The first canonical variate CV1 has discriminative power λ = 6.8485.The second CV2 has discriminative power equal to 0, however is helpful in perceiving the spread of the data points in the two groups.

Fig. 3 .

 3 Fig. 3. Kernel discriminant function obtained from Gaussian kernel with σ = 0.7. Negative and positive scores indicate the faulty and healthy state of the machine. Left panel: learning sample, first 300 faulty and next 300 healthy items. Right panel: testing sample, first 932 faulty and next 651 healthy items. Subsequent data items are linked by line segments. All items in both of the classes are classified perfectly.

gearboxes, the one in a faulty and the other in a healthy condition. The data were the basis of an elaboration presented in [START_REF] Bartelmus | A new feature for monitoring the condition of gearboxes in nonstationary operating systems[END_REF].

The raw recorded data (vibration series of the devices) were firstly segmented and transformed to frequency domain using power spectra densities. Next, 15 variables denoting amplitudes of derived spectra were obtained. This yielded two data matrices (A and B of size 1232 × 15, and 951 × 15 appropriately) characterizing the state of the faulty and the healthy device. The rows of the matrices constitute data vectors (called also data items) living in the 15-D space R 15 . Thus we have two samples of faulty and healthy data containing n1=1232 and n2=951 data items appropriately.

Taking as a new feature the sum of all the 15 variables, Bartelmus and Zimroz [START_REF] Bartelmus | A new feature for monitoring the condition of gearboxes in nonstationary operating systems[END_REF] were able to classify -on the base of the proposed feature -about 80 % of all data vectors. To classify the remainder, they needed an external variable indicating for the actual load of the working excavator.

The data were more thoroughly investigated in [START_REF] Zimroz | Investigation on spectral structure of gearbox vibration signals by principal component analysis for condition monitoring purposes[END_REF][START_REF] Bartkowiak | Outliers analysis and one class classification approach for planetary gearbox diagnosis[END_REF][START_REF] Zimroz | Multidimensional data analysis for condition monitoring: features selection and data classification[END_REF]. It appeared that the distribution of the variables is not Gaussian, the data contains a considerable number of outliers, moreover, the covariance structure in the two groups (faulty and healthy) is markedly different. It became evident, that a proper analysis of the data needs non-linear methods [START_REF] Zimroz | Multidimensional data analysis for condition monitoring: features selection and data classification[END_REF]. Now we have subdivided the entire data set into the learning and testing samples. The learning sample was obtained from randomly chosen 300 rows of the A matrix and other 300 rows of the B data matrix; this yielded the learning sample of size 600 × 15. The remainder of the data (932 rows from A and 651 rows from B) yielded the testing sample of size 1583 × 15.

Results

Results when using linear regression method

The linear regression was supposed to predict the class score given by the values: y i = -1 for data vectors i belonging to the faulty class, and y i = +1 for data vectors i belonging to the healthy class.

The calculations were done in Matlab. The estimates of the regression coefficients b = [ b1 , b2 , . . . , b15 , b0 ] T were obtained using the learning sample and the Matlab function mldivide, as described in Section 2.1.

The obtained estimates b were used to obtain the predicted values of the regression. This yielded the predicted vectors of discriminant scores ŷlearn for the learning and ŷtest for the testing samples. The obtained scores are visualized in Figure 1.

We have used the following decision rule: If for a data vector x i its predicted value satisfies ŷi ≤ 0, classify that data vector as 'faulty', otherwise -classify it as 'healthy'.

Looking at the results shown in Fig. 1 we find that: In the learning sample we have 3 incorrect assignments for a total of 600 items. The erroneous assignments happened in the faulty group, where 3 faulty items