$\DeclareMathOperator\argmax{argmax}$ $\DeclareMathOperator\argmin{argmin}$
Principal Component Analysis and K-means algorithm¶
PCA¶
Nous allons utiliser la PCA afin de réduire la dimension du jeu de données MNIST. Ceci va nous permettre de visualiser l'effet de la PCA. Pour ceux qui n'ont pas le jeu de données vous pouvez le téléchargez à l'adresse suivante : http://www.math.ens.fr/cours-apprentissage/mnist_digits.mat. Formez une matrice de design $X\in\mathbb{R}^{n\times d}$ en sélectionnant un sous ensemble de points (n=5000). Veillez toutefois à gardez un bon nombre de représentants par classe.
import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline
plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML, Markdown
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
'<script>'
'var waitForPlotly = setInterval( function() {'
'if( typeof(window.Plotly) !== "undefined" ){'
'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
'clearInterval(waitForPlotly);'
'}}, 250 );'
'</script>'
))
import urllib
# Interact widget
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from scipy.io import loadmat
from numpy.random import randint
mnist = loadmat('mnist_digits.mat') # données sous forme de dictionnaire
nbr_img = 5000 # nombre d'images à sélectionner au hasard
ind_rand = randint(len(mnist['x']), size=nbr_img) # indices des images sélectionnées
X, Y = mnist['x'][ind_rand, :], mnist['y'][ind_rand] # images, labels
# Histogramme des labels
trace = go.Histogram(
x=Y.flatten(),
xbins=dict(
start=-0.5,
end=9,
size=1
)
)
data = [trace]
layout = go.Layout(
title='Labels des images du jeu de données',
xaxis=dict(
title='Label'
),
yaxis=dict(
title='Nombre'
)
)
fig = go.Figure(data=data, layout=layout)
plotly.offline.iplot(fig)
1) Appliquez la PCA à la matrice de design $X$. Pour cela on pensera à tout d'abord centrer les données ($X_c=X-\frac{1}{n}\textbf{1}_nX^T\textbf{1}_n$). Ensuite on diagonalisera la matrice de corrélation : $X_c^TX_c$ (commande eig en Matlab).¶
Dans la suite on considérera que les vecteurs propres ont été ordonnés par valeurs propres décroissantes.
N, d = X.shape
n = 10**4
X_c = X - np.mean(X, axis=0)
C = X_c.T.dot(X_c)
eig, v = np.linalg.eigh(C)
eigenvalues = eig[-1], eig[-2]
vectors_PCA = v[:,-1], v[:, -2]
2) Tracez l'évolution du spectre de la matrice de corrélation. Rappelez l'interprétation de ces valeur propres et commentez ce que vous observez.¶
# Affichage de l'évolution de la distortion
trace_eig = go.Scatter(
x = list(range(d)),
y = eig[::-1],
mode = 'lines',
name = 'Valeurs propres'
)
layout_eig = go.Layout(
title= 'Évolution du spectre de la matrice de corrélation (échelle logarithmique)',
hovermode= 'closest',
xaxis= dict(
title= 'Itérations',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'Valeurs propres',
type='log',
ticklen= 5,
gridwidth= 2,
),
showlegend= True
)
plotly.offline.iplot(go.Figure(data=[trace_eig], layout=layout_eig))
3) Représentez sous forme d'image le premier et le second vecteur propre. Commentez.¶
plt.figure(figsize=(13,13))
plt.subplot(1, 2, 1)
plt.imshow(vectors_PCA[0].reshape(28, 28), cmap='gist_gray')
plt.title('Premier vecteur propre (associé au rayon spectral)')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(vectors_PCA[1].reshape(28, 28), cmap='gist_gray')
plt.title('Second vecteur propre')
plt.axis('off')
plt.show()
4) Projetez les images initiales sur l'espace engendré par les deux premiers vecteurs propres. Représentez alors le jeu de données MNIST sous forme d'un nuage de points.¶
Affichez différentes couleurs pour les différentes classes. En faisant le lien avec la question précédente, commentez la position des différents chiffres dans cet espace.
W = np.vstack(vectors_PCA).T
proj = X.dot(W)
# Affichage des données auxquelles K-means a été appliqué
traces_PDA = []
for i in range(10):
traces_PDA.append(go.Scatter(
x = proj[:, 0][Y.flatten() == i],
y = proj[:, 1][Y.flatten() == i],
mode = 'markers',
marker = dict(
size = 7,
showscale = False,
line = dict(
width = 1,
color = 'rgb(0, 0, 0)'
)
),
name = '{}'.format(i),
showlegend= True
))
layout_PDA = dict(
title = 'Projection des images sur l\'espace engendré par les deux premiers vecteurs propres',
xaxis = dict(
title = 'Premier vecteur propre',
ticklen = 5,
zeroline = False,
gridwidth = 2,
),
yaxis = dict(
title = 'Second vecteur propre',
ticklen = 5,
gridwidth = 2,
),
legend = dict(
orientation = 'h',
y = -0.2
)
)
plotly.offline.iplot(go.Figure(data=traces_PDA, layout=layout_PDA))
Pour y voir plus clair, affichons $0$ et $7$:¶
plotly.offline.iplot(go.Figure(data=[traces_PDA[0], traces_PDA[7]], layout=layout_PDA))
La composante des images représentant un $0$ :
- sur le premier vecteur propre, est globalement supérieure à celle des $7$. Cela s'explique, quand on regarde la visualisation de ce vecteur propre (cf. question précédente): vu comme une image, on distingue un zéro.
- sur le deuxième vecteur propre, est globalement inférieure à celle des $7$. De manière analogue, cela s'explique par le fait qu'on y distingue une forme ressemblant à un sept.
Un raisonnement similaire explique le positionnement relatif des autres chiffres dans l'espace de visualisation.
K-means¶
Commençons par quelques rappels sur l'algorithme. Soient $(x_i)_{i=1,\dots,n}$ des vecteurs appartenant à $\mathbb{R}^p$. On notera $X\in\mathbb{R}^{n\times p}$ la matrice de design associée. Le but est de partitionner ces données en $K$ clusters. On introduit alors une variable indicatrice $Z\in \lbrace 0,1\rbrace ^{n\times K}$ qui indique l'appartenance des points $x_i$ aux différentes clusters $k$. Plus précisément $Z_{ik}=1$ si $x_i$ est assignée au cluster $k$ et $Z_{ik}=0$ sinon. On introduit aussi $K$ ``moyennes" $(\mu_k)_{k=1,\dots,K}$ dans $\mathbb{R}^p$ (on verra dans la suite que les $\mu_k$ sont exactement les moyennes de chaque cluster). On notera $\mu\in\mathbb{R}^{K \times p}$ la matrice associée. Le critère qu'on utilise pour qualifier un bon partitionnement selon $Z$ et $\mu$ est la distorsion $J(\mu,Z)$ définie comme suit :
$$ J(\mu,Z) =\sum_{i=1}^n \sum_{k=1}^K Z_{ik}\Vert x_i-\mu_k \Vert^2. $$
Le but est alors de minimiser cette distorsion. L'algorithme K-means consiste en une minimisation alternée entre $\mu$ et $Z$ :
- Choix d'une matrice $\mu$ de manière aléatoire
- Etape 1: Minimisation de $J$ par rapport à $Z$ (revient à assigner chaque point $x_i$ au cluster le plus proche)
- Etape 2: Minimisation de $J$ par rapport à $\mu$ : $\mu_k = \frac{\sum_i Z_{ik}x_i}{\sum_i Z_{ik}}$
- Etape 3: Revenir à l'étape 1 jusqu'à convergence.
1) Montrer que l'étape 2 de l'algorithme correspond à l'annulation du gradient de $J$ par rapport à chacun des $\mu_k$ avec $Z$ fixé.¶
$$\frac{\partial J}{\partial μ_k}(μ, K) = \frac{\partial}{\partial μ_k} \left( \sum_{i=1}^n Z_{ik}\Vert x_i-\mu_k \Vert^2\right) = - 2 \sum_{i=1}^n Z_{ik} (x_i - μ_k)$$
Donc
$$ \begin{align*} & \quad \frac{\partial J}{\partial μ_k}(μ, K) = 0 \\ ⟺ & \quad \sum_{i=1}^n Z_{ik} x_i - μ_k \left(\sum_{i=1}^n Z_{ik}\right) = 0 \\ ⟺ & \quad \mu_k = \frac{\sum_i Z_{ik}x_i}{\sum_i Z_{ik}} \\ \end{align*} $$
2) Implémentez l'algorithme K-means (la fonction sqdist peut être à nouveau utile ici).¶
Celui-ci doit prendre en entrée une matrice de design X et un nombre de clusters K et sortir des moyennes $\mu$ ainsi que la matrice de correspondance Z.
Une bonne chose serait aussi de sortir l'évolution de la distortion au cours des itérations.
Remarque: Un critère d'arrêt de l'algorithme peut être lorsque la distortion n'évolue plus.
def K_means(X, K, returnClusterIndices=False):
n, p = X.shape
# Intialisation aléatoire de mu:
# Pour s'assurer que les moyennes sont distinctes,
# on les tire dans X où on a éliminé les lignes répétées
X_unique = np.unique(X, axis=0)
range_n = np.arange(len(X_unique))
np.random.shuffle(range_n)
mu = X_unique[range_n[:K]]
# liste de d'évolution de la distortion
distortions = []
while True:
# matrice des distances (au carré) deux-à-deux entre les x_i et les mu_i
# pairwise_dist[i, k] = distance entre x_i et mu_k
pairwise_dist = np.linalg.norm(X.reshape((n, 1, p)) - mu, axis=2)
# Etape 1: indices des clusters auxquels les x_i appartiennent
cluster_ind = np.argmin(pairwise_dist, axis=1)
# Calcul de la distortion courant et ajout à `distortions`
distortions.append(np.sum(pairwise_dist[np.arange(n), cluster_ind]))
# Etape 2: calcul des nouveaux mu_k
mu = np.array([np.mean(X[cluster_ind == k], axis=0) for k in range(K)])
# Condition d'arrêt: si la distortion n'évolue plus
if len(distortions)>1 and distortions[-1] == distortions[-2]:
break
Z = np.zeros((n, K))
Z[np.arange(n), cluster_ind] = 1
if returnClusterIndices:
return mu, Z, distortions, cluster_ind
else:
return mu, Z, distortions
3) Reprenez un des jeu de données du DM (par exemple accessibles à l'adresse www.di.ens.fr/~slacoste/teaching/apprentissage-fall2014/TP/classificationA.train).¶
Appliquez alors votre algorithme pour différentes valeurs de $K$. Tracez l'évolution de la distortion en fonction des itérations. Affichez ensuite les résultats (sous forme de nuage de points colorés en fonction du cluster).
letters = ['A', 'B', 'C']
values_of_K = {'A': 2, 'B': 2, 'C': 3}
ranges = {}
data = {}
trace = {}
mus = {}
layout = {}
for letter in letters:
ranges[letter] = {}
X_train_test = []
# Récupération des données d'entraînement et de test sous forme d'array
for train_test in ['train', 'test']:
with open('classification{}{}'.format(letter, train_test)) as f:
lines = []
for line in f:
lines.append(list(map(float, line.split())))
lines = np.array(lines).reshape([-1, 3])
X_train_test.append(lines[:, :2])
data[letter] = np.concatenate(X_train_test)
X_1 = data[letter][:,0]
X_2 = data[letter][:,1]
ranges[letter]['x_1'] = [int(np.min(X_1)-1), int(np.max(X_1)+1)]
ranges[letter]['x_2'] = [int(np.min(X_2)-2), int(np.max(X_2)+2)]
# K-Means
mu, Z, dist, cluster_ind = K_means(data[letter], values_of_K[letter], returnClusterIndices=True)
# Affichage des données auxquelles K-means a été appliqué
trace[letter] = go.Scatter(
x = X_1,
y = X_2,
mode = 'markers',
marker = dict(
size = 10,
color = cluster_ind,
showscale = False,
line = dict(
width = 2,
color = 'rgb(0, 0, 0)'
)
),
showlegend = False
)
layout[letter] = dict(
title = 'Ensemble de données {}'.format(letter),
xaxis = dict(
title = '$x_1$',
range = ranges[letter]['x_1'],
ticklen = 5,
zeroline = False,
gridwidth = 2,
),
yaxis = dict(
title = '$x_2$',
range = ranges[letter]['x_2'],
ticklen = 5,
gridwidth = 2,
),
legend = dict(
orientation = 'h',
y = -0.2
)
)
# Afficher les moyennes
mus[letter] = go.Scatter(
x = mu[:, 0],
y = mu[:, 1],
mode = 'markers',
marker = dict(
symbol = "star-diamond",
size = 15,
color = 'yellow',
line = dict(
width = 2,
color = 'black'
)
),
name = 'Moyennes',
showlegend= True
)
# Affichage de l'évolution de la distortion
trace_dist = go.Scatter(
x = list(range(1, len(dist)+1)),
y = dist,
mode = 'lines',
name = 'Distortion'
)
layout_dist = go.Layout(
title= 'Ensemble {}: Évolution de la distortion (échelle logarithmique)'.format(letter),
hovermode= 'closest',
xaxis= dict(
title= 'Itérations',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'Distortion',
type='log',
ticklen= 5,
gridwidth= 2,
),
showlegend= True
)
plotly.offline.iplot(go.Figure(data=[trace_dist], layout=layout_dist))
plotly.offline.iplot(go.Figure(data=[trace[letter], mus[letter]], layout=layout[letter]))
4) Nous allons maintenant regarder une application plus visuelle.¶
Choisissez une image de votre choix. On va alors utiliser K-means pour encoder l'image sur seulement $K$ couleurs. Pour cela on va considérer que chaque pixel $i$ de l'image est une observation $x_i$ appartenant à $\mathbb{R}^3$. Appliquez alors Kmeans en commençant avec $K=5$.
# Récupération de l'image
img = plt.imread(urllib.request.urlopen("http://younesse.net/images/app_stat/K-means/conan.png"))
# Pour se débarrasser de la transparence
img = img.dot(np.eye(4,3))
plt.imshow(img)
plt.show()
import time
start = time.time()
mu, Z, dist, cluster_ind = K_means(img.reshape([-1, 3]), 5, returnClusterIndices=True)
end = time.time()
print('Temps d\'exécution: {:.2f}s'.format(end - start))
# Attribution à chaque point de la couleur de la moyenne du cluster associé
cluster_ind_list = cluster_ind.tolist()
for ind, cluster in enumerate(cluster_ind_list):
cluster_ind_list[ind] = mu[cluster][:3]
img_new_colors = np.reshape(np.array(cluster_ind_list), img.shape)
plt.imshow(img_new_colors)
5) Affichez la nouvelle image ainsi obtenue (image qui contient pour chaque pixel la couleur associée à son cluster).¶
Essayez différentes valeurs de $K$.
def K_means_on_image(K, img=img):
mu, Z, dist, cluster_ind = K_means(img.reshape([-1, 3]), K, returnClusterIndices=True)
# Attribution à chaque point de la couleur de la moyenne du cluster associé
cluster_ind_list = cluster_ind.tolist()
for ind, cluster in enumerate(cluster_ind_list):
cluster_ind_list[ind] = mu[cluster][:3]
return np.reshape(np.array(cluster_ind_list), img.shape)
display(Markdown("### K-Means appliqué à l'image pour différentes valeurs de K"))
nrows, ncols = 3, 3
plt.figure(figsize=(13,13))
for k in range(2, ncols * nrows+2):
plt.subplot(nrows, ncols, k-1)
plt.imshow(K_means_on_image(k))
plt.title('K = {}'.format(k))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
@interact(K=widgets.IntSlider(min=1,max=20, step=1, value=3))
def show_K_means_on_image(K):
plt.figure(figsize = (7., 7.))
plt.imshow(K_means_on_image(K))
plt.title('K = {}'.format(K))
plt.axis('off')
plt.show()
display(Markdown("### Version interactive"))
6) Bonus (pour aller plus loin) :¶
Comme vous l'avez peut être remarqué l'algorithme K-means dépend beaucoup de son initialisation et il peut s'avérer assez lent. Afin de régler ce problème, Arthur et al. ont proposé une meilleure initialisation (qui possède en plus des garanties théoriques quant au vrai écart à l'optimum). Cette méthode a été baptisé K-means ++ et consiste à remplacer l'étape 0 de l'algorithme K-means par l'algorithme suivant :
- Etape 0: Choisir de manière uniforme un point $x_i$ comme étant le premier centre.
- Etape 1: Pour chaque point $x_i$ calculez la distance entre ce point et le centre le plus proche. On note $D_{\mu_t}(x_i)$ cette distance.
- Etape 2: Choisir un nouveau point $x_i$ comme étant un nouveau centre mais cette fois ci en suivant une loi de probabilité pondéré par $D_{\mu_t}(x_i)^2$.
- Etape 3: Répétez l'étape 1 et 2 jusqu'à avoir $K$ centres.
Implémentez cette méthode et vérifiez qu'elle fonctionne mieux en pratique (meilleur objectif final ainsi que convergence plus rapide).