Chi square test Cross Table: Example (2/2) Total per Category Cat A1 (30) Cat A2 (30) Cat A3 (40) Cat B1 (40) 12 12 16 Cat B2 (30) 9 9 12 Cat B3 (30) 9 9 12

Categorical variables : Chi-squared test Let us consider two variables i and j so that i ∈ [1..r ] and j ∈ [1..c]. Let o i,j be the number of observed data so that the first feature's value is i and the second feature's value is j. Let n i be the total number of elements having i as a value for there first feature. Let N be the total number of observations. Then, the expected contingency e i,j computes as follows:

e i,j = n i • n j N Computing the χ 2 χ 2 = r i=1 c j=1 (o i,j -e i,j ) 2 e i,j
If there is no link, the repartition in the cross table should be almost exactly proportional with the size of the categories.

If the observed values o i,j are far from this supposed proportional repartition, there is probably a link between some of the categories.

Cross Table : Example (3/3) The Chi-squared is a test of independence. Therefore, the hypothesis that you are trying to reject is H 0 : "The two variables are independent".

Computing the p-value

The p-value is computed using the known distribution of the Chi-squared function (using tables, a calculator, or a graph) considering degrees of freedom of the problem.

Degrees of freedom = (r -1)(c -1)

A p-value close to zero rejects the hypothesis that the two variables are independent.

We can say that there is "1-p-value % chance" that the link between the two variables is not random luck.

Computing further indexes will be necessary to assess the degree of correlation.

One way to start explaining a chi2 result is to look at which associations have a strong influence on the chi2.

This can be done by analyzing the normalized residuals:

r ij = o ij -e ij √ e ij
And their contribution to the chi2, which is very useful for large cross-tables:

c ij = 100 × r 2 ij χ 2 J. Sublime Data Analysis -Lecture 3 ISEP 2021/2022 13 / 36 It is denoted T ∈ [0, 1].
It measures the amount of dependency between two categorical variables.

Chuprov coefficient

T = χ 2 N (c -1)(r -1)
If T is close to 0, we can't predict one from the other.

If T is close to 1, the link is strong and we can predict one from the other.

Cramer's V is another is another index measuring the amount of dependency between two variables.

Cramer's V V = χ 2 N • min(c -1, r -1)
If V is close to 0, we can't predict one from the other.

If V is close to 1, the link is strong and we can predict one from the other.

While the p-value evaluates whether a result based on the chi-squared has good chances of being significant, its value is not proportional to the amount of correlation.

Chuprov coefficient is best reliable with square contingency tables, while Cramer's V works best with rectangular ones.

Both Chuprov coefficient and Cramer's V can be very biased: unevenly distributed observations, one variable with much more possible values than the other, etc.

Both Chuprov coefficients and Cramer's V can't be used if the result of the p-value shows that the chi-squared result is not significant.

The Chi-squared test is not recommended with very small data sets (most expected values bellow 10), and can be replaced by the similar Fisher test in such cases.

A significant dependency between two variables using the chi squared test does not imply causality. The chi squared result alone with the p-value are not informative enough:

Chuprov or Cramer's V are needed to assess strength of the link and determine the reliability of predicting one variable based on the other. Computing the residuals and contribution to the chi2 explains how the different pairing of categories affect the test and might thus be remarkable or related.

Example: Dataset

T = χ 2 N (c -1)(r -1) = 10.67 95 √ 2 = 0.28 V = χ 2 N • min(c -1, r -1) = 10.67 95 • min(1, 2) = 0.34

Interpretation

From the p-value, we deduce that the two variables show a significant dependency.

On a rectangle contingency matrix Cramer's V indicates around 34% of correlation between the two variables (so does the Chuprov coefficient with 28%).

We have seen that for numerical variables we can use the covariance, the correlation (Spearman or Pearson), and the determination coefficient.

For categorial variables, we have the Chi square (or Fischer) test coupled with a Cramer's V or Chuprov coefficient.

What about mixed data ?

What are mixed data ? Mixed data are datasets that contains both numerical and categorial data (and sometimes other things).

Numerical bivariate analysis won't work on them (because of the categorial variables)

Neither should the chi square test (because of the continuous variables)

From numerical to categorial data We see that a few of the observed values that are too high or too low compared with the expected values.

Example

Example: Contribution to the chi 2

Let us take now a look at the contribution to the chi2: In addition to the previous informations, we clearly see that haplogroups I and R1a have the most influence on the chi2, which means that these 2 groups have an influence on the size:

c ij A-H I J L-Q R1a Other
It explains why the chi2 is significant But it also explain why the link remains weak since only 2 haplogroups are remarkable, which far from enough to predict the size from the haplogroup, or the opposite.

  Given this contingency table, is there a correlation between the Y haplogroup and the average height ?

					Example Example Example Example Example
	From numerical to categorial data Example: Dataset Example: Groups of heights Example: Summing Example: Chi 2 Example: Chi 2 Example: residual analysis
		A-H	I	J	L-Q	R1a	Other
	Dealing with mixed data ? 175 Dealing with mixed data ? 176	41 64	13 14	7 5	26 17	0 0	40 60
	Numerical bivariate analysis won't work on them (because of the Numerical bivariate analysis won't work on them (because of the 177 64 15 0 21 1 61 178 64 13 11 23 8 58 To better explain this result, we can take a look at the normalised
	categorial variables) Neither should the chi square test (because of the continuous variables) categorial variables) Neither should the chi square test (because of the continuous variables) So, how do we deal with them ?! Discretizing the numerical variables 200 0 3 0 0 0 4 Numerical variables can be categorized by building artificial groups 179 69 7 10 22 7 70 180 53 15 5 25 11 54 181 61 0 6 17 3 60 182 62 20 3 18 8 63 183 52 17 2 13 11 50 184 38 21 2 11 9 39 185 53 20 3 10 18 52 186 47 11 0 6 12 46 187 38 17 2 11 0 39 188 33 0 0 0 3 3 189 0 7 0 10 9 20 190 0 7 1 0 7 10 193 28 0 1 5 0 0 195 20 8 0 4 5 13 196 7 9 0 4 0 11 197 6 0 1 0 1 0 A-H I J L-Q R1a Other 175-179 302 62 33 109 16 289 180-184 266 73 18 84 42 266 185-189 171 55 5 37 42 160 190-194 28 7 2 5 7 10 195-199 33 17 1 8 6 24 200+ 3 5 0 2 2 A-H I J L-Q R1a Other e i,j A-H I J L-Q R1a Other e i,j A-H I J L-Q R1a Other residuals: n j 175-179 302 62 33 109 16 289 811 180-184 266 73 18 84 42 266 749 185-189 171 55 5 37 42 160 470 190-194 28 7 2 5 7 10 59 195-199 33 17 1 8 6 24 89 200+ 3 5 0 2 2 6 18 175-179 296.55 80.88 21.79 90.48 42.47 278.83 180-184 273.88 74.70 20.12 83.56 39.22 257.51 185-189 171.86 46.87 12.63 52.44 24.61 161.59 190-194 21.57 5.88 1.59 6.58 3.09 20.28 195-199 32.54 8.88 2.39 9.93 4.66 30.60 200+ 6.58 1.80 0.48 2.01 0.94 6.19 175-179 296.55 80.88 21.79 90.48 42.47 278.83 r ij A-H I J L-Q R1a Other 180-184 273.88 74.70 20.12 83.56 39.22 257.51 175-179 0.32 -2.10 2.40 1.95 -4.06 0.61 185-189 171.86 46.87 12.63 52.44 24.61 161.59 190-194 21.57 5.88 1.59 6.58 3.09 180-184 -0.48 -0.20 -0.47 0.05 0.44 0.53 20.28 195-199 32.54 8.88 2.39 9.93 4.66 185-189 -0.07 1.19 -2.15 -2.13 3.50 -0.13 30.60 200+ 6.58 1.80 0.48 2.01 0.94 190-194 1.38 0.46 0.33 -0.62 2.22 -2.28 6.19 195-199 0.08 2.73 -0.90 -0.61 0.62 -1.19 6 n i 803 219 59 245 115 755 2196 d f = (6 -1) × (6 -1) = 25 χ 2 = 87.33, 200+ -1.40 2.39 -0.70 -0.01 1.09 -0.08
	202 between ranges of interest 3 203 0	0 0	0 0	2 0	0 1	1 1
	207 The main difficulty is to find the right size/number of 0 1 0 0 1 0 212 0 1 0 0 0 0
	groups/categories
	J. Sublime	Data Analysis -Lecture 3	ISEP 2021/2022	31 / 36

  Let us take now a look at the contribution to the chi2:

	c ij	A-H	I	J	L-Q	R1a	Other
	175-179 0.11 5.05 6.61 4.34 18.89 175-179 0.11 5.05 6.61 4.34 18.89	0.42 0.42
	180-184 0.04 0.26 0.00 180-184 0.26 0.04 0.26 0.00	0.23 0.23	0.32 0.32
	185-189 0.00 1.61 5.28 5.20 14.07 185-189 0.00 1.61 5.28 5.20 14.07	0.02 0.02
	190-194 2.19 0.24 0.12 0.44 190-194 2.19 0.24 0.12 0.44	5.67 5.67	5.97 5.97
	195-199 0.01 8.52 0.93 0.43 195-199 0.01 8.52 0.93 0.43	0.44 0.44	1.63 1.63
	200+ 200+	2.23 6.55 0.55 0.00 2.23 6.55 0.55 0.00	1.36 1.36	0.01 0.01

Example

Example: Contribution to the chi 2

J. SublimeData Analysis -Lecture 3 ISEP 2021/2022

Example: Chi 2 interpretation

With χ 2 = 87.83 and d f = 25, we have a p-value bellow 0.01: There is a significant link between the Y haplogroup and the height of an individual ! We also have T = V = 0.089, so the link appears to be weak.

J. Sublime

Data Analysis -Lecture 3 ISEP 2021/2022

Example: chi squared table

With χ 2 = 10.67 and d f = 2, we have a p-value bellow 0.01: There is a significant correlation between hair color and eye color. What about mixed variables ?

We have seen that for numerical variables we can use the covariance, the correlation (Spearman or Pearson), and the determination coefficient.

For


Correlation & Regression

Keywords:

Regression and time series analysis: to visualize the data and guess the model Clustering: To guess the number of clusters and their shape, or to detect anomalies/outliers. Classification: To assess classes' separability.

Why visualizing data ?

A better visualization leads to a better understanding.

Graphical visualization is an efficient way to communicate information clearly and efficiently.

A lot of data are high dimensional and difficult to grasp without a lower dimensionality projection. Many Machine Learning techniques can benefit from data visualization:

Examples

Univariate data

Represents the distribution of a variable

Represents categorical variable distribution

Represents trends and error margins The curse of dimensionality

The "curse of dimensionality " is a term introduced in 1961 by Bellman refering to the problem of the explosive increase in data volume associated with adding extra dimensions in a mathematical space.

Why the curse of dimensionality matters

High dimensional data are difficult to efficiently visualize in 2D without loosing meaningful information.

While computer algorithms can grasp information in higher dimensional spaces than humans, past some point, they too have trouble extracting the meaningful information.

We will illustrate this problem with a simple example.

The Toy problem

We consider a Pattern Recognition problem with 3 classes. We have 9 available observation represented in 1 dimension.

A simple approach would be to do the following:

Divide the feature space into uniform bins.

Compute the class ratio in each bin and use a majority vote to classify the bin.

When a new example comes, put it in its closest bin.

The Toy problem 2D

If we add a 2nd dimension we pass from 3 cases in 1D to 3 2 = 9 in 2D.

We have an another problem: do we maintain the density of examples per bin or do we keep the same number of example than in 1D?

The Toy problem 2D

If we add a 2nd dimension we pass from 3 cases in 1D to 3 2 = 9 in 2D.

We have an another problem: do we maintain the density of examples per bin or do we keep the same number of example than in 1D?

If we want to maintain the density, we now need 27 examples instead of 9.

Choosing to maintain the number of examples to 9 will result in a very sparse 2D scatter plot.

The Toy problem 3D

Adding a 3rd dimension makes the problem worse:

We move to 3 3 = 27 bins.

Keeping a constant density now requires 81 observations.

Keeping the same number of observation gives us an almost empty 3D plot.

The curse of dimensionality

The curse of dimensionality generates several phenomenons such as:

The concentration of observations in given space areas

The desertification of the data space

The depopulation of the center of hyper-volumes

In practice, the curse of dimensionality means that, for a given sample size, there is a maximum number of variables beyond which a classifier performance will degrade rather than improve.

The approach used on the toy example is mostly ineffective:

There are other approaches less affected by the curse of dimensionality, but the problem will still exist and will always re-surface at some point. Exponential growth of the complexity of the target function (which estimates the density). To learn well, the target function requires a dense enough learning space.

For one dimension, we have plenty of different density functions. However, for high dimension, we only have the Gaussian multivariate density. For large values of D, the density can only be treated in Gaussian simplified forms (due to the inverse covariance matrix in the Gaussian mixture model).

How to deal with the curse of dimensionality ?

These findings suggest that we need special treatment to manipulate high dimensional data:

Incorporating prior knowledge to reduce the search space.

Providing smoother and more generic target functions.

Reducing dimensionality.

Reducing dimensionality

The reduction of dimensionality by the mean of space transformations, projections, or feature selection is usually the prefered solution in combination with the two others.

Prior knowledge is not always available, and may alone result in over-fitting.

Smoother target functions alone may result in a global decrease of performances and are cumbersome to design.

Dimension Reduction

Dimension reduction is based on the hypothesis that high dimensional data are not uniformly distributed:

There are high density and low density areas.

Dimension reduction looks for structures that define these high and low density areas to tackle the curse of dimensionality.

Dimension Reduction

Dimension reduction is based on the hypothesis that high dimensional data are not uniformly distributed:

There are high density and low density areas.

Dimension reduction looks for structures that define these high and low density areas to tackle the curse of dimensionality.

There are many techniques of dimension reduction:

Linear and non-linear

Deterministic and probabilistic

Supervised and unsupervised

Selection vs Extraction

There are two main methodology for dimension reduction:

Feature selection: Choosing a subset of all the features.

Feature extraction: Creating a small set of new features by combinations of the original ones.

        X 1 X 2 . . . X D         feature selection ---------→     X i 1 X i 2 X i M             X 1 X 2 . . . X D         feature extraction ----------→     Y 1 Y 2 Y M     = f                 X 1 X 2 . . . X D                 Outline Introduction Feature selection methods
Linear Feature extraction methods Bibliography Q1 How to measure the relevance of the variables ?

We need to find a measure of relevance, or evaluation criterion J(X ), to quantify the importance of one variable or a combination of variables.

Feature selection

Q1 How to measure the relevance of the variables ?

We need to find a measure of relevance, or evaluation criterion J(X ), to quantify the importance of one variable or a combination of variables.

Q2 How to obtain the optimal subset ?

To answer this question, we need to define a procedure, or algorithm, that will determine the optimal subset of relevant variables based on the evaluation criterion.

Feature selection

Q1 How to measure the relevance of the variables ?

We need to find a measure of relevance, or evaluation criterion J(X ), to quantify the importance of one variable or a combination of variables.

Q2 How to obtain the optimal subset ?

To answer this question, we need to define a procedure, or algorithm, that will determine the optimal subset of relevant variables based on the evaluation criterion.

Q3 Which optimality criterion to use ?

With the answer to Q2 being an optimization algorithm, we need to define a stopping criterion for this optimization process.

Choosing a criterion for feature selection: Examples

For a Classification problem, we can test the discriminant quality of the system in the presence presence or absence of a variable.

For a Regression problem, we can test the quality of the prediction with respect to the other variables.

Feature selection: Branch & Bound based methods

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Feature selection: Branch & Bound based methods

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Feature selection: Branch & Bound based methods

Branch & Bound based methods use an optimization strategy that finds the optimal subset of variable based on a minimal bound to reach on the evaluation criterion.

Branch & Bound methods can only work if the selection criterion J(X ) is monotonic.

Monotonic criterion: Definition Let us consider two subsets of features A 1 and A 2 . Then,

J(•) is monotonic if and only if A 1 ⊂ A 2 =⇒ J(A 1 ) ≤ J(A 2 ).
Branch & Bound methods can only work if the selection criterion J(X ) is monotonic.

Monotonic criterion: Definition Let us consider two subsets of features A 1 and A 2 . Then,

J(•) is monotonic if and only if A 1 ⊂ A 2 =⇒ J(A 1 ) ≤ J(A 2 ).
Problem: Most effective evaluation criteria are not monotonic ... The use of sub-optimal methods is then required:

Sequential Forward Selection (SFS) Sequential Backward Selection (SBS)
Bidirectional Selection (BS)

Sequential Forward Selection

Let X = {X 1 , • • • , X D } be a set of variables.
The SFS procedure is the following:

Initially set the selected set of variables A 0 as empty.

At each step k, select the variable X i that maximizes the following criterion of evaluation J(A k ):

J(A k ) = max X i ∈(X \A k-1 ) J(A k-1 ∪ X i )
This will result in a list of variables ordered by their importance. The procedure can be stopped at any step when a given criterion (e.g. acceptable loss of information) has reached an acceptable value.

Sequential Backward Selection

Let X = {X 1 , • • • , X D } be a set of variables. The SBS procedure is the following:

Initially start with a set A D = X containing all the variables.

At each step for

k ∈ [D -1 • • • 1],
remove the variable X i of the least importance according to:

J(A k ) = max X i ∈A k+1 J(A k+1 \X i )
This will also result in a list of variables ordered by their importance, with the most important being at the end of the list. The procedure can also be stopped has soon as the remaining set is below an acceptable lower bound of accuracy on a given criterion.

Bidirectional Selection

The Bidirectional Selection procedure performs the search in both Forward and Backward direction in a competitive manner:

The procedure stops when one of the following case occurs:

One of the two direction has found the best subset of variables before reaching the middle of the search space.

The two branches are reaching the middle.

This method reduces the search time as the search is performed in both directions and stops when there is a solution regardless of the direction.

Remark

The sets of best selected variables found respectively by SFS and SBS are not equal and rarely identical because of their different principles of selection.

Introduction

Feature selection methods

Linear Feature extraction methods

Bibliography

Introduction

Feature Extraction: Principle

Feature extraction consists in building new features from the original ones with one or several of the following goals:

Having a lower number of features while keeping a maximum of information

Having better features with which the data are easier to process

Having feature with which the data are easier to visualize

       X 1 X 2 . . . X D        feature extraction ----------→     Y 1 Y 2 Y M     = f               X 1 X 2 . . . X D              
Reducing the dimensionality: goals

The Principal Component Analysis (PCA) is a method of data analysis and dimensionality reduction that seeks the directions of space that best represent the correlations between D random variables.

Goals

1 Obtaining an approximation of a scatter of K objects in a subspace of low dimension (2 or 3).

2 Summarizing a data set represented by a matrix X of N rows and D columns.

3 Getting a simple and reliable visualization of the information contained in the data.

Reducing the dimensionality: Principle of PCA

The Principal Component Analysis (PCA) is a method of data analysis and dimensionality reduction that seeks the directions of space that best represent the correlations between D random variables.

Principle 1 Compute the covariances or correlations between variables.

2 Use this indexes to create new variables:

The original variables are replaced with the new variables.

Each new variable is a linear combination of the original variables. These new variables must be as independent (uncorrelated) from each other as possible. PCA and Normalizing the data: Centering and reducing

PCA : Principle

Prior to a Principal Component Analysis, centering the variables around their means and scaling them may be necessary.

Centered variables

M =    x 1,1 -µ 1 • • • x 1,D -µ D . . . . . . . . . x N,1 -µ 1 • • • x N,D -µ D    Centered Reduced variables M =     x 1,1 -µ 1 σ 1 • • • x 1,D -µ D σ D . . . . . . . . . x N,1 -µ 1 ) σ 1 • • • x N,D -µ D σ D    
PCA and the Variance-Covariance Matrix

Once the data have been normalized, computing the variance-covariance matrix is the next step in order to assess the potential relationships between the different variables.

Variance-Covariance

C =    VAR(X 1 ) • • • COV (X 1 , X D ) . . . . . . . . . COV (X D , X 1 ) • • • VAR(X D )   

Calculation of new PCA variables

Now that we know the covariance relationships between all variables, we seek to create new variables using linear combinations of the old ones:

A first new variable represented by a single axis u 1 so that the projection of u 1 has a maximum variance.

Then we seek for a second axis u 2 independent of u 1 (i.e. orthogonal to u 1 ) that best explains the remaining variance.

And so on.

Calculation of new PCA variables

Now that we know the covariance relationships between all variables, we seek to create new variables using linear combinations of the old ones:

A first new variable represented by a single axis u 1 so that the projection of u 1 has a maximum variance.

Then we seek for a second axis u 2 independent of u 1 (i.e. orthogonal to u 1 ) that best explains the remaining variance.

And so on.

Note that for visualization purposes, 2 or max 3 axes are usually enough.

PCA : Calculating the new variables (2/6)

We seek u 1 of norm 1 so that the projection of the data on u 1 has a maximum variance.

Computing u 1

Let us note V π the variance of π u (M):

V π = 1 N -1 π u (M) T • π u (M) = u T • M T M N -1 C •u = u T • C • u
With C the variance-covariance matrix of the original dataset X .

Note that it only works because all variables M i are centered ! PCA : Calculating the new variables (3/6)

Computing u 1 From V π = u T • C • u, maximizing V π is equivalent to solve: C • u = V π • u
By definition, the solution of this equation is the set of couples : eigenvectors / eigenvalues of C .

The largest eigenvalue λ 1 is the variance V π , the maximum reachable.

The associated eigenvector is the projection axis u 1 along which V π is maximum.

PCA : Calculating the new variables (4/6)

The above reasoning has allowed us to conclude that the vector which explains the more inertia (i.e. the variance) of the data points is the first eigenvector. Similarly, the second vector which explains most of the remaining inertia is the second eigenvector, etc.

We also saw that the variance of the data projected onto the n th eigenvector (also called -variance explained by this axis) is λ n .

Remark

In Matrix form, we have: C ∝ UΛU T where Λ is a diagonal matrix containing all single values.

PCA : Calculating the new variables (5/6)

Computing u k

To find the k th component, the data are updated:

Xk = X - k-1 s=1 X • u s • u T s
The process that we have just presented is repeated by finding the new maximal eigenvectors of the variance-covariance matrix of the new data.

It turns out that this process gives the remaining eigenvectors of the original variance-covariance matrix.

Therefore, in practice the k th component is picked by choosing the k th largest eigenvalue and the associated eigenvector.

Principal Component Analysis

PCA : Calculating the new variables

Another possible algorithm

1 Choose P the proportion of information (also called inertia) that we want to keep.

2 Select y eigenvectors so that:

y = argmin y y i=1 λ i λ ≥ P
We know that the new variables represent linear combinations of the old ones. It is therefore interesting to analyze the correlation between the old and new variables to interpret the results.

Correlations between an old variable and a new variable

r (M i , π j ) = 1 N T M i • π j σ i λ j
This expression follows directly the original formula for the correlation coefficient, recalling that λ j is the variance of the projection of the data on the new u j axis.

Circle of correlations (1/3)

Now that we know the correlation between the old and new variables, we can simply visualize a circle of correlations.

This allows to both characterize the new variables (which makes it possible to visually interpret the data) and to visualize the correlation between the old variables.

Principle SVD Given a centered data matrix M, the goal of the SVD is to find the matrices U, S and V to achieve the following decomposition: With the SVD, we have M = USV T , therefore:

C = VSU T USV T N -1 = V S 2 N -1 V T
Thus:

The eigenvalues can be computed from the singular values λ i =

s 2 i N-1
V stores the eigenvectors of M T M.

Besides dimension reduction, MDS is well suited for applications or algorithms based on the distances between the data, rather than the data themselves:

Re-building the map of a country based on the distance between the cities.

Building a 2D map from 3D coordinates (e.g. star constellation map).

etc.

Remarks

With N(N-1) 2 distances, it is always possible to generate the position of N points in a N dimension space.

MDS computes an approximation for a lower dimensional space.

We have

The points

x 1 , • • • , x N in D dimensions
Any distance δ ij between points x i and x j

We want to find The points y 1 , • • • , y N in 2 or 3 dimensions

The distance d ij between any y i and y j so that it is close to

δ ij
Given that d ij is a function of y i and y j (the δ ij are constant), thereafter are a few possible objective functions:

J aa = i<j (d ij -δ ij ) 2 i<j δ 2 ij
penalizes large absolute errors

J rr = i<j d ij -δ ij δ ij 2
penalizes large relative errors

J ar = 1 i<j δ ij i<j (d ij -δ ij ) 2 δ ij A compromise between the two
The last one, J ar , is also referred as Sammon Criterion.

Figure :

 : Figure: From 21 variables to 2 variables using PCA. Note that the 3 clusters are still visible.

MM

  The data centered data matrix of size N × D U The matrix containing the left singular vectors (eigenassay), of size N × D. S The diagonal matrix containing the singular values, of size D × D V T The matrix containing the right singular vectors (eigengenes), of size D × D Relation to PCA Let C be the variance-covariance matrix of M. Then: C = M T M N -1

  Exponential growth of the number of examples required to maintain a given sample density: For a density of N examples/bin in D dimensions, N D examples are required.

J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

The algorithm starts with all features and removes one at each step.J. SublimeData Analysis -Lecture 4 ISEP

2021/2022

The algorithm starts with all features and removes one at each step.

If the evaluation criterion scores below the bound, the algorithm backtracks and tries to remove another feature.J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

The algorithm starts with all features and removes one at each step.

If the evaluation criterion scores below the bound, the algorithm backtracks and tries to remove another feature.

Groups of features below the bound are remembered so that sub-sets of these groups are never tested later in the process.J. SublimeData Analysis -Lecture

ISEP 2021/2022

The algorithm starts with all features and removes one at each step.

If the evaluation criterion scores below the bound, the algorithm backtracks and tries to remove another feature.

Groups of features below the bound are remembered so that sub-sets of these groups are never tested later in the process.

When there is no feature left to try removing, the algorithm remembers the best score and backtracks to another branch.J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

The algorithm starts with all features and removes one at each step.

If the evaluation criterion scores below the bound, the algorithm backtracks and tries to remove another feature.

Groups of features below the bound are remembered so that sub-sets of these groups are never tested later in the process.

When there is no feature left to try removing, the algorithm remembers the best score and backtracks to another branch.

The algorithm stops when all branches with valid subsets have been tested.J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

Select among the new variables the ones that best represent the distribution of the data and eliminate the others.J. SublimeData Analysis -Lecture

ISEP 2021/2022

Data Analysis -Lecture 4 ISEP 2021/2022

Choose p the number of dimensions that you want to keep.

Select the eigenvectors (u 1 , • • • , u p ) associated with the p largest eigenvalues.

Project the data on the new axes in order to obtain a description of these data according to the new variables (called principal components): π u (M) = M • u.

If the number of dimensions is ≤ 3, it is possible to represent the data in the new space.J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

Project the data on the new axes.

If the number of dimensions is ≤ 3, it is possible to represent the data in the new space.J. SublimeData Analysis -Lecture 4 ISEP 2021/2022

The Toy problem

We consider a Pattern Recognition problem with 3 classes. We have 9 available observation represented in 1 dimension.

Feature selection

Definition

Variable selection is a process to choose an optimal subset of relevant variables from a set of variables, according to a performance criterion.

We can ask three basic questions: Q1 How to measure the relevance of the variables ? Q2 How to obtain the optimal subset ? Q3 Which optimality criterion to use ? We seek u 1 of norm 1 so that the projection of the data on u 1 has a maximum variance.

Computing u 1

Principal Component Analysis PCA : Calculating the new variables (2/6)

We seek u 1 of norm 1 so that the projection of the data on u 1 has a maximum variance.

Computing u 1

Let us note V π the variance of π u (M):

Principal Component Analysis PCA : Calculating the new variables (2/6)

We seek u 1 of norm 1 so that the projection of the data on u 1 has a maximum variance.

Computing u 1

Let us note V π the variance of π u (M):

With C the variance-covariance matrix of the original dataset X .

J. Sublime Data Analysis -Lecture 4 ISEP 2021/2022

PCA : Calculating the new variables (5/6)

Computing u k

To find the k th component, the data are updated:

The process that we have just presented is repeated by finding the new maximal eigenvectors of the variance-covariance matrix of the new data.

J. Sublime Data Analysis -Lecture 4 ISEP 2021/2022

PCA : choosing the number of components

We typically decide on the number of principal components required by examining a scree plot. Look for a point at which the proportion of variance explained by each subsequent principal component drops off, this is often referred to as an elbow in the scree plot.

J. Sublime Data Analysis -Lecture 4 ISEP 2021/2022

Reminders

Reminder: eigenvalues The eigenvalues of a matrix X are determined by the vector λ so that:

Where | • | is the determinant operator.

Reminder: eigenvectors

The eigenvector v i linked to a eigenvalue λ i solves the following equation:

Limits of PCA

The problem with PCA Using PCA for large data sets with many attributes ?

Examples: genetic data, text data The covariance matrix is of size d × d. If there are too many attributes, computing the eigenvalues is too computationaly expensive because the matrix is too big.

Solution: Singular Value Decomposition (SVD)

Powerful algorithms exist in Python, R, matlab and SAS to compute it.

When to use SVD ?

More attributes than data: N < D When the centered data matrix M contains mostly 0 (sparse data). Multi-dimensional scaling: Cost function

We must search d ij so that it minimizes an objective function.

We can define such function in a general manner:

Multi-dimensional scaling: Example 2

The dataset "Eurodist" represents the distance in kilometers between 21 European cities.

The source data set is a square matrix containing the distances between the cities:

J. Sublime Data Analysis -Lecture 4 ISEP 2021/2022

Multi-dimensional scaling: Conclusions

The different MDS algorithms differ in the following points:

The distance used in the source space.

The Stress (objective) functions: Using a different stress function will lead to a different result.

The optimization procedure: Linear MDS can be solved analytically but cannot model complex (non-linear) low dimensional manifolds well. However, non-linear MDS requires a more complex and iterative algorithm.


Dr Jérémie Sublime
email: jeremie.sublime@isep.fr

Data Analysis -Lecture 5 Data visualization: non-linear methods

Limits of linear methods (2/2) Linear methods such as PCA cannot discover the structure of a data set with a non-linear shape. Figure: Example of a spiral data set that could be modeled in 2D or even 1D, but that PCA and non-linear methods cannot handle.

Plan

ISOMAP: Principle

The Isometric feature mapping (Isomap) [Tenebaum et al. 2000] is a reduction method that can handle non-linear subspaces. It is based on the following principles:

For neighboring samples (close data points), the Euclidian distance provides a good approximation of the geodesic distance.

For distant points, the geodesic distance can be approximated using an Euclidian distance-based neighborhood graph of the data and searching the shortest path between the two points.

1 Build a distance neighborhood graph G using the Euclidian distance in the input space: This can be done using KNN, or by connecting the points within a radius .

2 Approximate the geodesic distance between all data points: Compute the shortest path between all points of the graph using Dijkstra algorithm or Floyd algorithm.

3 Apply the Multi-dimensional scaling (MDS) algorithm to the geodesic distance matrix from Step 2.

Remark

The Isomap algorithm can be seen as a non-linear variant of the MDS algorithm.

ISOMAP: Algorithmic complexity

The main issue of the Isomap algorithm is that it can be quite slow. For a data set of size N, an input space of dimension D, an output space of dimension d, and considering a neighborhood of size k in the first step, we have:

Step

1: K-nearest neighbors O(N 2 D)
Step

2: Djikstra O(N 2 log (N) + N 2 k), Floyd O(N 3 )
Step 3: MDS algorithm O(N 2 d)

Isomap

MDS reminders

Multi-dimensional scaling: Algorithm

Choose an objective function J and a step parameter 0 < η < 1.

Compute the distance matrix δ if it is not provided.

Initialize the projected points randomly and compute the corresponding projected distance matrix.

Optimize the position of the projected points until convergence using gradient descent

The "Swiss roll" Example (1/3)

The "swiss roll" data set contains 20000 points.

We will show thereafter the development of the Isomap algorithm on this data set. Locally Linear Embedding: Principle

The locally Linear Embedding algorithm (LLE) [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] addresses the same problem as Isomap in a different way:

It preserves the local properties of the data by representing each point by a linear combination of its nearest neighbors.

It builds a projection to a linear low dimensional space preserving the neighborhood.

Locally Linear Embedding

Algorithm

Locally Linear Embedding: Principle

The locally Linear Embedding algorithm (LLE) [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] addresses the same problem as Isomap in a different way:

It preserves the local properties of the data by representing each point by a linear combination of its nearest neighbors.

It builds a projection to a linear low dimensional space preserving the neighborhood.

LLE algorithm

1 Compute the k nearest neighbors

2 Compute the weights needed to reconstruct each point using a linear combination of its neighbors

3 Project the results in a lower dimensional space using the weights

Locally Linear Embedding: Steps Locally Linear Embedding: Computing the weights (1/2)

The local geometry is modeled by linear weights that reconstruct each data point as a linear combination of its nearest neighbors.

Let us note W = (w ij ) N×N the weight matrix, where w ij is the weight given to x j in the reconstruction of x i .

Reconstruction error cost function

(W ) = N i=1 |x i - j =i w ij x j | 2
The weights w ij are minimized subject to two constraints:

1 Each data point is reconstructed only from its neighbors (i.e. w ij = 0 is x i and x j are not neighbors)

2 Rows of the weight matrix W sum to 1:

∀i j w ij = 1
Locally Linear Embedding: Computing the weights (2/2)

Let us consider a sample x i with its k nearest neighbors η ij . With the constraint j w ij = 1, the weights can be found in 3 steps:

Locally Linear Embedding

Algorithm

Locally Linear Embedding: Computing the weights (2/2)

Let us consider a sample x i with its k nearest neighbors η ij . With the constraint j w ij = 1, the weights can be found in 3 steps:

1 Compute the neighborhood correlation matrices C jk and their inverse:

C jk = (x i -η ij ) • (x i -η ik )
Locally Linear Embedding: Computing the weights (2/2)

Let us consider a sample x i with its k nearest neighbors η ij . With the constraint j w ij = 1, the weights can be found in 3 steps:

1 Compute the neighborhood correlation matrices C jk and their inverse:

C jk = (x i -η ij ) • (x i -η ik )
2 Compute the lagrangian multiplier λ that enforces the constraint

j w ij = 1: λ = 1 -j,k C -1 jk (x T i η ik ) j,k C -1 jk
Locally Linear Embedding: Computing the weights (2/2)

Let us consider a sample x i with its k nearest neighbors η ij . With the constraint j w ij = 1, the weights can be found in 3 steps:

1 Compute the neighborhood correlation matrices C jk and their inverse:

C jk = (x i -η ij ) • (x i -η ik )
2 Compute the lagrangian multiplier λ that enforces the constraint

j w ij = 1: λ = 1 -j,k C -1 jk (x T i η ik ) j,k C -1 jk 3
Compute the reconstruction weights:

w ij = k C -1 jk (x T i η ik + λ)
Locally Linear Embedding: Mapping the data (1/2)

Once the weights W have been computed, we seek to find the

Y = {y 1 , • • • , y n }, y i ∈ R d the new points in a low dimensional space.
This can be done by searching Y that minimizes an embedding cost function similar to the reconstruction error cost function from the input space.

Embedding error cost function

φ(Y ) = N i=1 |y i - j =i w ij y j | 2
Locally Linear Embedding

Algorithm

Locally Linear Embedding: Mapping the data (2/2)

The embedding vectors Y i are found by minimizing the cost function φ(Y ).

To make the optimization problem well-posed we introduce two constraints: j y j = 0 and

1 N i Y i Y T i = I
Locally Linear Embedding

Algorithm

Locally Linear Embedding: Mapping the data (2/2)

The embedding vectors Y i are found by minimizing the cost function φ(Y ).

To make the optimization problem well-posed we introduce two constraints: j y j = 0 and

1 N i Y i Y T i = I Then, the cost function becomes: φ(Y ) = i,j M ij (Y T i Y j ) Where M ij = δ ij -W ij -W ji + k w ki w kj δ ij is equal to 1 if i = j and 0 otherwise
Locally Linear Embedding: Mapping the data (2/2)

The embedding vectors Y i are found by minimizing the cost function φ(Y ).

To make the optimization problem well-posed we introduce two constraints: j y j = 0 and

1 N i Y i Y T i = I Then, the cost function becomes: φ(Y ) = i,j M ij (Y T i Y j ) Where M ij = δ ij -W ij -W ji + k w ki w kj δ ij is equal to 1 if i = j and 0 otherwise
The optimal embedding for a d-dimensional space is found by computing the bottom d + 1 eigenvectors of the matrix M.

Locally Linear Embedding

Example

Locally Linear Embedding: Face data set Example

The initial data represent images of faces.

In the 2D space, these images are grouped according to the position, lighting and expression. It is a nonlinear dimensionality reduction technique well-suited for embedding high-dimensional data for visualization in a low-dimensional space of two or three dimensions.

Stochastic Neighbor Embedding (SNE)

SNE: Basic idea Encoding high dimensional neighborhood information as a distribution

Intuition: Random walk between data points with a high probability to jump to a close point.

Finding low dimensional points such that their neighborhood distribution is similar.

How to measure the distance between distribution ?

The most common measure to do so is the KL divergence SNE: local embedding algorithm (1/4)

Consider the neighborhood around any input data point

x i ∈ R D
Imagine that we have a Gaussian distribution centered around x i Then the probability that x i chooses some other datapoint x j as its neighbor is in proportion with the density under this Gaussian.

A point closer to x i will be more likely than one further away t-distributed stochastic neighbor embedding SNE algorithm SNE: local embedding algorithm (2/4) P j|i the probability that x i chooses x j as neighbor is computed as follows:

P j|i = exp(-||x i -x j || 2 /2σ 2 i ) k =i exp(-||x i -x k || 2 /2σ 2 i )
P i|i = 0 σ i the chosen size of the neighborhood:

a low σ i will lead to a smaller local neighborhood with stronger probabilities a high σ i leads to uniform weights σ i can be set up differently for each point.

The final distribution over pairs is symmetrized:

P ij = 1 2N (P j|i + P i|j ) t-distributed stochastic neighbor embedding SNE algorithm SNE: local embedding algorithm (3/4) Given a dataset X = {x 1 , • • • , x N } ∈ R D , we define the distribution P ij Goal: Finding a good embedding Y = {y 1 , • • • , y N } ∈ R d , d << D,
and ideally d should be 2 or 3.

Measuring and minimizing the embedding quality

For all point in Y we can define Q in a similar way than P (notice the absence of σ i and the assymetry):

Q ij = exp(-||y i -y j || 2 ) k =i exp(-||y i -y k || 2 )
We optimize Q to be close to P by minimizing the KL-divergence on the parameters

Y = {y 1 , • • • , y N } ! The KL divergence (1/2)
Measuring the distance between two discrete distributions P and Q:

KL(Q||P) = i Q(i) log Q(i) P(i)
It is a metric function as it is not symmetric. We have P and we are looking for

Y = {y 1 , • • • , y N } ∈ R d such that the distribution Q we infer will minimize L(Q) = KL(P||Q) (notice the Q on the right, which is uncommon). Note that KL(P||Q) = ij P ij log P ij Q ij = -ij P ij (Q ij ) + const We can show that ∂L ∂y i = j (P ij -Q ij )(y i -y j )
Remarks This is not a convex problem: No convergence guarantees ! Multiple restarts may be needed.

Another issue may arise: crowding problems.

t-distributed stochastic neighbor embedding SNE algorithm

Crowding problems

We have seen at the beginning of this class that in high dimension we have more room: it means that points can have a lot of different neighbors in a given radius.

However, in lower dimension (1D, 2D or 3D), there is a lot less room in the same radius. This is the "crowding problem": we might not have enough room to accommodate all neighbors. This is one of the biggest weakness for LLE and SNE.

t-SNE solution to crowding problems

Change the Gaussian in Q by a heavy tailed distribution : if Q changes slower, we have more "wiggle room" to place points.

(x) ≈ 1 + x 2 v -(v +1)/2 For v = 1, we get: p(x) ≈ 1 1+x 2
The probability of this distribution goes much more slower to zero than a Gaussian. This is convenient for our crowding problem.

We can show that this is equivalent to averaging a Gaussian with some prior over σ 2 . Then, we can redefine Q ij so that:

Q ij = (1 + ||y i -y j || 2 ) -1 k l =k (1 + ||y i -y k || 2 ) -1
t-distributed stochastic neighbor embedding t-SNE algorithm

t-SNE gradients

We can show that the gradients of t-SNE objective function are:

∂L ∂y i = j (P ij -Q ij )(y i -y j ) 1 + ||y i -y j || 2
Compared with the SNE gradients:

∂L ∂y i = j (P ij -Q ij )(y i -y j )
credit: "Visualizing Data using t-SNE"

Both repulse close dissimilar points and attract far similar points, but the t-SNE has a smaller attraction term that reduces crowding issues. We can define the desired perplexity and set σ i to get that (bisection method): values between 5 and 50 usually work well.

It is an important parameter as it allows to capture different scales in the data.

Neural networks

Neural networks are a variety of algorithms that mimic brain cells:

Each neuron, or each group of neuron specializes in recognizing certain patterns. This incoming information is not only received by a single neural cell, but also influences other cells in its neighborhood.

Self-organizing maps (SOM) [T. Kohonen 1984] are a family of unsupervised neural networks based on this principle:

It learns in an unsupervised way using prototypes that play the role of neurons.

The prototypes follow a given topology so that neighborhood dependencies exist between them.

Groups of neighbor prototypes specialize in recognizing similar groups of data. Self-organizing maps: Architecture (2/2)

C (W ) = N i=1 K k=1 K k,χ(x i ) ||x i -w k || 2
χ(x i ) is the id of the best matching prototype w j so that: χ(x i ) = argmin j ||x i -w j || 2 . 0 ≤ K k,j ≤ 1 is the neighborhood topological function (e.g. K k,j = 0 if w j and w k are far on the topological grid).

Self-organizing maps: Topological learning Remark Isomap and LLE try to find a low dimensional manifolds in the data using neighborhood dependencies of the data.

SOM has the opposite approach by trying to wrap a low-dimensional manifold around the data using prototypes neighborhood dependencies.

Self-organizing maps for visualization purposes

Once they are built, SOM can be used in several ways for visualization purposes:

Cluster visualization on the 2D map using the best matching neuron Topological visualization of the data set using colors to project distance matrices between the prototypes on the map.

Component plane analysis by projecting the components on the map.

Data projection and transformation

Self-organizing maps: Visualization Example 3 The initial topology of the map does not matter much: If the map is big enough, it will work.

Limitations of SOM

The learning rate is a critical parameter.

There are plenty of variations of the algorithm each of them with different strengths and weaknesses.

Introduction Isomap

Locally Linear Embedding t-distributed stochastic neighbor embedding Self-organizing maps Bibliography

  Euclidian distance VS Geodesic distanceEuclidian based distances are one of the main problem when dealing with data embedded in a non-linear subspace.Custom distances embedded in the low-dimensional manifold, such as the geodesic distance, can help dealing with such data.

Figure

  Figure: Euclidian distance Figure: Geodesic distance

Figure

  Figure: A 1000 points extract of the Swiss roll data set, with an example of Euclidian and Geodesic distances.

Figure :

 : Figure: LLE algorithm [Roweis and Saul 2000]

Figure :

 : Figure: Images placed at the bottom of the figure correspond to successive points encountered on the line at the top right, sweeping a continuum of facial expression. [Roweis and Saul 2000]

  .02596 -0.12176 + 0.24408 = 0.09637 The KL divergence (2/2) KL properties KL(Q||P) ≥ 0 and is a convex function KL(Q||P) = 0 if and only if Q = P if P(i) = 0 but Q(i) > 0, then KL(Q||P) = ∞ credit: https://timvieira.github.io/blog/post/2014/10/06/kl-divergence-as-an-objective-function/ t-distributed stochastic neighbor embedding SNE algorithm SNE: local embedding algorithm (4/4)

Figure :

 : Figure: Example of a single layer neural network

Figure :

 : Figure: Example of a rectangle grid SOM architecture with 4 neighbors per neuron: Each example activates a winning neuron and also affects its neighbors.

Figure :

 : Figure: Number recognition using SOM [Rogovschi et al. 2008]

  We have seen that σ i is a key neighborhood parameter. For each distribution P j|i (depending on σ i ), we define the perplexity:Perp(P j|i ) = 2 H(P j|i ) where H(P) = -

	t-distributed stochastic neighbor embedding	t-SNE algorithm
	t-SNE: Perplexity (1/2)
		P i log(P i )
		i
	If P is uniform over k elements -perplexity is k.
	Smooth version of k in kNN
	Low perplexity means small σ 2
	High perplexity means large σ 2

t-SNE algorithm source : "Visualizing Data using t-SNE"

J. SublimeData Analysis -Lecture 5 ISEP 2020/2021

Data Analysis -Lecture 5 ISEP 2020/2021

ISEP 2020ISEP /2021

Self-organizing maps: Algorithm (1/2)

Initialization step

Define the topology of the map Randomly initialize the prototypes for each neuron Define a learning rate function (t)

Competition step

Choose a data x i randomly Determine the winning neuron χ(x i ) = argmin 1≤j≤K ||x i -w j || 2 Adaptation step: Update the prototypes

Repeat 2 and 3 until the prototypes updates are insignificant.

Self-organizing maps: Visualization Example 1


Data partitioning : Examples We want to automatically discover groups of similar data

How many clusters ?

Guessing the right number of clusters is one of the oldest problem in clustering, and several methods require it as a parameter.

Clusters, similarity and distance

The similarity or dissimilarity between two observations is most often computed using a distance function.

Depending on the data and the application, some distances may prove more appropriate than others (e.g. Hamming distance for text data).

Euclidian distance ||a -b|| 2 = i (a i -b i ) 2 Squared Euclidian distance ||a -b|| 2 2 = i (a i -b i ) 2 Manhattan distance ||a -b|| 1 = i |a i -b i | Maximum distance ||a -b||∞ = max i |a i -b i | Mahalanobis distance (a -b) S -1 (a -b)
where S is the covariance matrix Hamming distance

Hamming (a, b) = i (1 -δ a i ,b i )

Table: Examples of common distances

Creating custom distance functions is sometimes required.

Clustering partitions

Let us denote X = {x 1 , • • • , x N }, x n ∈ R d a data data set containing N observations described by d real features, and C = {c 1 , ..., c K } the set of possible clusters.

Clustering partitions

The result of a clustering algorithm is called a clustering partition, or just partition. Depending on whether or not a data can belong to several clusters, the partition can be hard, soft, or fuzzy.

(a) Hard clustering

c 1 c 2 c 3 x 1 1 0 0 x 2 0 1 0 x 3 0 0 1 x 4 0 0 1 (b) Soft clustering c 1 c 2 c 3 x 1 1 1 0 x 2 0 1 1 x 3 0 0 1 x 4 0 0 1 (c) Fuzzy clustering c 1 c 2 c 3 x 1 0.9 0.1 0 x 2 0 0.8 0.2 x 3 0 0.3 0.7 x 4 0 0 1.0
c 1 c 2 c 3 x 1 1 0 0 x 2 0 1 0 x 3 0 0 1 x 4 0 0 1 (b) Soft clustering c 1 c 2 c 3 x 1 1 1 0 x 2 0 1 1 x 3 0 0 1 x 4 0 0 1 (c) Fuzzy clustering c 1 c 2 c 3 x 1 0.9 0.1 0 x 2 0 0.8 0.2 x 3 0 0.3 0.7 x 4 0 0 1.0
Let us note S the clustering partition:

In hard clustering, S is a vector

S = {s 1 , • • • , s N }, s n ∈ [1..K ] In soft clustering, S is a matrix S = {s n,k } (N×K ) , s n,k ∈ (0, 1)
In fuzzy clustering, S is a matrix too:

S = {s n,k } (N×K ) , s n,k ∈ [0, 1] ∀n, k s n,k = 1

Types of clustering algorithms

Now that we know what we want to do, the question is: Given a data set, a distance function and a desired type of partition, how do we find the clusters ?

There are several "families" of clustering algorithms that favor different methods to extract the clusters from the data:

Density

Spectral clustering algorithms

Sees clustering as a graph partitioning problem using the similarity matrix between the different objects.

Examples of algorithms: Shi-Malik algorithm (normalized cuts), Meila-Shi algorithms.

Mainly use for image analysis

Properties Pros: Good results, does not assume any shape for the clusters, fast with sparse data.

Cons: High complexity O(N 3 ), does not scale well, not very intuitive. The "Density" is based on the number of data within the radius "Eps" of a given element.

Outline

A point is a core point if it has more data in its neighboring radius than a specified number of points "MinPts".

Core points are the main structures for the clusters.

A border point has fewer points in its radius Eps than the number MinPts but it is within range of a core point.

A noise point has fewer points in its radius Eps than the number MinPts and is not within range of a core point.

DBSCAN: Limits

Difficult to parametrize High Computational cost Does not handle well varying densities OPTICS is a good evolution of DBSCAN that handles varying densities a lot better.

HCA: Introduction

Principle Create a partition a each step by aggregating together the 2 closest clusters.

Each data starts out as a single cluster.

During the process a cluster can be a single data, or a group of data.

The algorithm return a hierarchy of partitions in the form of a tree containing the history of the different aggregations.

HCA: Algorithm

1 Repeat until there is only one cluster left:

1.a Compute the distance matrix between all existing clusters 1.b Merge the two closest clusters 1.c Update the dendrogram (hierarchical tree)

2 Cut the tree according to a criterion of your choice to get the final partition.

This algorithm requires distance criterion between clusters: the linkage criterion.

HCA: Linkage criteria

Name Formula Comments

Single-linkage

Ds (c 1 , c 2 ) = min x∈c 1 ,y ∈c 2 d(x, y ) Distance between the two closest elements Complete-linkage Dc (c 1 , c 2 ) = max x∈c 1 ,y ∈c 2 d(x, y ) Distance between the two farthest element Average-linkage Da(c 1 , c 2 ) = 1 |c 1 ||c 2 | x∈c 1 y ∈c 2 d(x, y ) Average pairwise distance Centroid-linkage Dµ(c 1 , c 2 ) = ||µ 1 -µ 2 ||

Centroid-linkage

A good compromise between single and complete link.

Less sensitive to noise and outliers.

tends to favor elliptical clusters.

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022

HCA: Other linkages

Average-linkage

Less noise sensitive, tends to favor hyper-spherical clusters High computational cost for roughly the same results as the centroid-linkage.

Centroid-linkage

A good compromise between single and complete link.

Less sensitive to noise and outliers.

tends to favor elliptical clusters.

Ward's method

This similarity is based on the increase on the mean squared distance to the centroid when 2 clusters are merged.

Less sensitive to noise and outliers.

tends to favor elliptical clusters.

K-Means: Introduction

The method of the K Means, or K-Means algorithm, is a special case of the mobile centers methods. Its principle consists in using a certain number of representatives (centroids, or prototypes) in the data space. Each of these representative will represent a cluster.

In the end, each data is associated to its closest prototype in order to obtain the clustering partition.

The groups formed with this process are homogeneous and well separated.

1 Randomly select K initial center.

2 Assign each data x to its closest center µ i .

3 If the partition doesn't change: stop 4 Else, update the centers, each µ i must be the gravity center of its cluster:

µ i = 1 |c i | x∈c i x 5 Go to 2
Like in most clustering algorithms, step 2 is dependent on a chosen distance function. For the K-Means algorithm, the Euclidian distance is usually prefered:

d(x, y ) = d i=1 (x i -y i ) 2
Guaranteed to monotonically decrease average squared distance in each iteration

L(µ) (t) = N i=1 min j ||µ j -x i || 2 2 L(µ) (t+1) ≤ L(µ) (t)
Convergence to a local minimum Algorithmic complexity at each iteration :

O(n • d • K )
Guaranteed to monotonically decrease average squared distance in each iteration This problem can be reduced by running the algorithm several time.

L(µ) (t) = N i=1 min j ||µ j -x i || 2 2 L(µ) (t+1) ≤ L(µ) (t)
It is also possible to initialize the centers using another clustering algorithm.

EM algorithm for the GMM: introduction

The Expectation-Maximization algorithm (EM) is an iterative method used to find the maximum likelihood or the maximum a posteriori estimate of parameters in probabilistic and statistical models.

Adapatation to clustering

In clustering, it is possible to represent clusters using parametric distribution: Each cluster has its own distribution with its own parameters.

The EM algorithm can then be used to figure out these parameters and link each data to a cluster.

EM algorithm for the GMM: Gaussian mixture model

Gaussian model in dimension 1 µ ∈ R the mean value of the gaussian distribution σ ∈ R + the standard deviation of the distribution

N (µ, σ) = 1 σ √ 2π exp -1 2 (x-µ) 2 σ 2
Gaussian model in dimension d µ ∈ R d the mean value of the gaussian distribution Σ = (σ i,j ) d×d the variance-covariance matrix of the distribution

N (µ, Σ) = exp -1
EM algorithm for the GMM: Optimization process

The goal of the EM algorithm is to find the parameters Θ that maximize the likelihood of generating the observed data:

Θ ML = Argmax Θ log P(X |Θ) = Argmax Θ N i=1 P(x i |Θ) Θ ML = Argmax Θ N i=1 log K k=1 π K • N (µ k , Σ k , x i )
This is done via an iterative two-step process: Expectation step: Update the latent variables (cluster labels) using the current parameters Θ: assign each data to the cluster c that maximizes p(x ∈ c).

Maximisation step: Use the new latent variables to update the parameters: Update the mean value, covariance matrix and mixing probabilities of all clusters.

EM algorithm for the GMM: Update rules

E-Step P(x n ∈ c) = s n,c = π c N (µ c , Σ c , x) K k=1 π k N (µ k , Σ k , x) = 1 Z π c exp -1 2 (x-µc ) T Σ -1 c (x-µc ) |Σ c |(2π) d M-Step N c = N n=1 s n,c =⇒ π c = N c N µ c = 1 N c N n=1 s n,c • x n Σ c = 1 N c N n=1 s n,c • (x n -µ c )(x n -µ c ) T
EM algorithm for the GMM: Algorithm

Algorithm 1: EM Algorithm for the GMM Initialize Θ randomly while the partition S is not stable do Picking a similarity measure Whatever the type of clustering algorithm, they all rely on a similarity measure. Picking a good similarity measure is therefore of paramount importance.

E-Step: evaluate S = argmax S p(S|X , Θ) forall x n ∈ X do s n (k) = 1 Z π c N (x n , µ c , Σ c ) end M-Step: Re-evaluate Θ forall c ∈ [1..K ] do Update µ c , Σ c
It is the core of any model.

It determines what type of structures are considered interesting.

It will ultimately change the clustering result.

Remarks

Choosing a similarity measure, or building one often require a good knowledge of the data to analyze, or at least of the field they come from.

Picking a similarity measure already introduces a bias in a task that is supposed to be exploratory.

Picking a similarity measure

Validating clustering results

Validating clustering results is a difficult process:

There is no proper definition of what a "good cluster" is, or looks like.

There is no "ground truth" in unsupervised learning to check the quality of a partition.

Clusters come in all shapes and forms.

The number of clusters to be found is usually unknown, and sometimes there are none.

Not all data sets have well-defined clusters.

Identifying clusters

While it is usually obvious to visually spot the clusters in 2D or 3D data sets, it is almost impossible to do so when there are more features !

Internal indexes

Internal indexes are criteria that can be used to evaluate the quality of a clustering partition.

Internal indexes can be subjective as they favor some cluster shapes, and can be biased toward a lower number of clusters.

They usually assess the compactness of the clusters and whether they are well separated.

Remark

Internal indexes are said to be "internal", by opposition to "external indexes" that compare a result with an external ground truth.

Internal indexes: Davies-Bouldin Index Let S i be the average scatter of a cluster c i around its mean value µ i :

S i = 1 |c i | x∈c i ||x -µ i || 2
Let M i,j = ||µ i -µ j || 2 be the average distance between two clusters c i and c j .

Davies-Bouldin Index for K clusters

DB = 1 K K i=1 max j =i S i + S j M i,j
Internal indexes: Davies-Bouldin Index

The Davies-Bouldin index has the following properties:

A lower DB value means a better clustering.

This index is not normalized.

It favors spherical clusters.

it is biased so that it gives lower values with less clusters.

Remark

The Davies-Bouldin index is usually the favored internal index used to evaluate partitions created using the K-Means algorithm.

Internal indexes: Silhouette Index

The Silhouette index is another internal index with interesting properties: It is normalized between -1 and 1 (a value bellow 0 meaning a bad partition).

It can evaluate if a a data belongs to a cluster, if a cluster is well formed, or the whole partition. Let a x ||x -µ i || 2 be the mean distance between x ∈ c i and the data that belong to the same cluster. And b x the mean distance to the data that belong to the closest other clusters.

b x = min k =i 1 |c k | y ∈c k ||x -y || 2 Silhouette index of an element x SC (x ∈ c i ) = b x -a x max(a x , b x )

Internal indexes: Calinski & Harabasz

This index is defined as follows: Criteria to pick the right number of clusters Some criteria exist to pick the optimal number of clusters, but they require to run the algorithm several times with different number of clusters in order to evaluate the solutions:

CH(k) = B/(k -1) W /(N -k) B is
The Bayes Information Criterion (BIC) and Akaike Information Criterion (AIC) work well with probabilistic algorithms.

The Minimum Description Length criterion (MDL) and the Minimum Message Length criterion (MML) are recommended for hierarchical clustering.

Stability can also be used to guess the right number of clusters.

Remark

Most of the time, the best solution is to ask the advice of an expert in the field of the studied data.

Picking a clustering algorithm

Picking the right clustering algorithm depends on several factors and requires prior knowledge on the data. With high dimensional data sets, you usually don't have this knowledge:

Are there any clusters ?

What is the shape of the clusters and how well separated are they ?

Other considerations may include the complexity of your algorithm:

Complexity depending on the size of the data set.

Complexity depending on the number of attributes of the data set.

Picking an algorithm often ends up being a trade-off between the expected quality of the result and the computational complexity of your algorithm.

Figure :

 : Figure: Example of 2 clusters with Gaussian distributions

  not presume the shape of the clusters or their number, detect outliers. Cons: High complexity O(N 2 ), difficult to parametrize Hierarchical clustering algorithms Inspired from phylogenetic trees in biology. Approaches: agglomerative, divisive Examples of algorithms: HCA, CURE, CLINK Properties Pros: Highlight hierarchical cluster structures, comprehensive visualization of the clusters Cons: High complexity O(N 2 log N)toO(N 3 ), choosing where to cut the dendrogram can be cumbersome. Prototype-based clustering algorithms Uses vector quantization to sum up the data using representatives (e.g. mean values) Examples of algorithms: K-Means, Fuzzy C-Means Properties Pros: Lower complexity O(N) Cons: Requires to give the number of clusters and to presume of their shapes. Distribution-based clustering algorithms Uses strong mathematical models to describe the clusters (e.g. mixture of gaussian distributions) Examples of algorithms: EM Properties Pros: Lower complexity O(N), strong models, convergence proofs Cons: Requires to give the number of clusters and to assume that the clusters will follow a given distribution.

  Less noise sensitive, tends to favor hyper-spherical clusters High computational cost for roughly the same results as the centroid-linkage.

  Convergence to a local minimum Algorithmic complexity at each iteration : O(n • d • K ) Non-convex optimization : NP-hard problem K-Means: Limits Instability The K-Means algorithm is highly dependent on the initialization of the centers.

Figure :

 : Figure: Depending on the similarity measure, both partitions can make sense.

  the sum of square distances between clusters W is the sum of square distances within clusters k is the number of clusters Properties of the CH-Index Not normalized Better when higher With balanced clusters, the CH index is generally a good criterion to indicate the correct number of clusters.

J. SublimeData Analysis -Lecture 6 ISEP 2021/2022

(x-µ) T Σ -1 (x-µ) |Σ|(2π) d J. SublimeData Analysis -Lecture 6 ISEP 2021/2022

DBSCAN: Algorithm

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022

K-Means: Properties

Guaranteed to monotonically decrease average squared distance in each iteration

Convergence to a local minimum Clusters with different density or sizes can also be difficult to find using the K-Means algorithm.

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022

K-Means: Limits

Choosing the number of clusters The K-means algorithms requires the number of clusters "K " to be provided as a parameter. In purely exploratory data mining tasks, this number is not always known and must be guessed.

Sometimes the only solution is to try several values for K and to see which gives the best results.

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022

EM algorithm for the GMM: Gaussian mixture model

The Gaussian mixture model consists in modeling clusters as Gaussians described by their mean value, their covariance matrix and a mixing probability.

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022

EM algorithm for the GMM: Gaussian mixture model

The Gaussian mixture model consists in modeling clusters as Gaussians described by their mean value, their covariance matrix and a mixing probability.

On note P(x i |Θ) la probabilité d'observer x i sachant les paramètres de tous les clusters, avec EM algorithm for the GMM: Limits

The EM algorithm for the GMM model detects only ellipsoid clusters.

Computing Σ -1 the inverse of the variance-covariance matrix can be both difficult and time consuming in high dimension.

Can be solved by using diagonal matrices, but it is less accurate.

Very much like the K-Means algorithm, the K needs to be provided and the EM algorithm is dependent on its initialization.

Can be improved by initializing the EM algorithm using K-Means or HCA ... EM algorithm for the GMM: Limits

The EM algorithm for the GMM model detects only ellipsoid clusters.

Computing Σ -1 the inverse of the variance-covariance matrix can be both difficult and time consuming in high dimension.

Can be solved by using diagonal matrices, but it is less accurate.

Very much like the K-Means algorithm, the K needs to be provided and the EM algorithm is dependent on its initialization.

Can be improved by initializing the EM algorithm using K-Means or HCA ...

Remark

The K-Means algorithm is a degenerate case of EM algorithm for the GMM where Σ = I d and ∀c, π c = 1 K .

Internal indexes: Silhouette Index

Silhouette index of a partition

The silhouette index best value is 1.

It favors spherical clusters. Using External indexes to asses cluster quality Quite often clustering algorithms will be tested in a non-exploratory setting:

Range

A data set for which the real classes are known will be used for test purposes (with the added advantage that the number of clusters will be known)

In this case, the cluster found will be compared to the real classes using external indexes or purity measures.

Cluster purity

Examples of external indexes

Rand Index, Adjusted Rand Index, Accuracy, etc.

Summary of open problems in clustering

Defining the notion of similarity, and picking or defining a distance function accordingly.

Picking the right algorithm and the right parameters for this algorithm.

Figuring out how many clusters to search for.

Evaluating clustering results.

Comparing clustering results.

All known data analysis problems: normalizing the data or not, picking relevant attributes, removing the outliers or not, etc.

J. Sublime Data Analysis -Lecture 6 ISEP 2021/2022