Versions¶
In [1]:
#version de Python
import sys
sys.version
Out[1]:
'3.12.7 | packaged by Anaconda, Inc. | (main, Oct 4 2024, 13:17:27) [MSC v.1929 64 bit (AMD64)]'
In [2]:
#version de numpy
import numpy as np
np.__version__
Out[2]:
'1.26.4'
In [3]:
#numba
import joblib
joblib.__version__
Out[3]:
'1.4.2'
Chargement et préparation des données¶
In [4]:
import sklearn.datasets as ds
covtype = ds.fetch_covtype()
In [5]:
#matrice des descripteurs
X = np.array(covtype.data)
X.shape
Out[5]:
(581012, 54)
In [6]:
#attribut classe
target = np.array(covtype.target)
target.shape
Out[6]:
(581012,)
In [7]:
#valeurs possibles
np.unique(target,return_counts=True)
Out[7]:
(array([1, 2, 3, 4, 5, 6, 7]),
array([211840, 283301, 35754, 2747, 9493, 17367, 20510],
dtype=int64))
In [8]:
#standardisation des X
from sklearn.preprocessing import StandardScaler
std = StandardScaler()
Z = std.fit_transform(X)
#dim.
Z.shape
Out[8]:
(581012, 54)
Régression logistique - Fisher Scoring¶
In [88]:
# code proposé par COPILOT, à la question :
#
# Peux-tu me proposer un code Python pour le calcul de la régression
# logistique binaire avec l'algorithme Fisher scoring ?
"""
def fit_reg(Z, y, max_iter=5, tol=1e-6):
# copie locale
X = np.copy(Z)
# ajouter la constante
#X = np.hstack([np.ones((X.shape[0], 1)), X])
X = np.concatenate((np.ones((X.shape[0],1)),X),axis=1)
# Initialisation des poids
beta = np.zeros(X.shape[1])
for i in range(max_iter):
# Calcul des probabilités prédictives
p = 1.0/(1.0 + np.exp(-1.0 * X @ beta))
# Calcul de la matrice de la pente (gradient)
gradient = X.T @ (y - p)
# Calcul de la matrice Hessienne - argh la taille de W
#W = np.diag(p * (1 - p))
#H = X.T @ W @ X
#solution moins gourmande en espace mémoire
proba = p * (1 - p)
M = X.T * proba
H = M @ X
# Mise à jour des poids avec Fisher scoring
beta_new = beta + np.linalg.inv(H) @ gradient
# Condition d'arrêt
if np.sum(np.abs(beta_new - beta)) < tol:
break
beta = beta_new
return beta
"""
Régresion logistique binaire¶
In [9]:
#régression logistique binaire
#max_iter est passé à 10 par rapport à la vidéo "Numba"
def fit_reg(Z, y, max_iter=10, tol=1e-6):
# copie locale
X = np.copy(Z)
# ajouter la constante
X = np.concatenate((np.ones((X.shape[0],1)),X),axis=1)
# Initialisation des poids
beta = np.zeros(X.shape[1])
for i in range(max_iter):
# Calcul des probabilités prédictives
p = 1.0/(1.0 + np.exp(-1.0 * X.dot(beta)))
# Calcul de la matrice de la pente (gradient)
gradient = X.transpose().dot(y - p)
#calcul de la matricie Hessienne
#solution moins gourmande en espace mémoire
proba = p * (1 - p)
M = X.T * proba
H = M.dot(X)
# sécurisation - rajouter une constante dans la diagonale
for j in range(H.shape[0]):
H[j,j] += 1e-6
# Mise à jour des poids avec Fisher scoring
beta_new = beta + np.linalg.inv(H).dot(gradient)
# Condition d'arrêt
if np.mean(np.abs(beta_new - beta)) < tol:
break
beta = beta_new
return beta
Régression one-vs-rest -- Non-parallélisée¶
In [10]:
#on part du principe que les classes sont identifiées par 1, 2, 3, ....
def one_vs_rest(Z, cible, max_iter=10, tol=1e-6):
nb_classes = np.max(cible)
#pour chaque modalité de la variable cible
lst_coefs = []
for k in range(nb_classes):
y = (cible==(k+1))
lst_coefs.append(fit_reg(Z,y,max_iter,tol))
#return liste des coefficients
return lst_coefs
In [11]:
#entraînement sur nos données
coefs_ = one_vs_rest(Z,target)
print(len(coefs_))
7
In [12]:
#coefs sous-forme matricielle
betas_ = np.array(coefs_)
betas_.shape
Out[12]:
(7, 55)
Fonction de prédiction¶
In [13]:
# fonction de prédiction
def pred_reg(Z,betas):
# copie locale
X = np.copy(Z)
# ajouter la constante
X = np.concatenate((np.ones((X.shape[0],1)),X),axis=1)
# produit matriciel
logit = X @ betas.T
#max. par ligne
id_max = np.argmax(logit,axis=1)
# prediction
return id_max + 1
In [14]:
# appliquer
pred = pred_reg(Z,betas_)
np.unique(pred,return_counts=True)
Out[14]:
(array([1, 2, 3, 4, 5, 6, 7], dtype=int64),
array([205641, 305579, 49147, 1222, 438, 3711, 15274],
dtype=int64))
In [15]:
# accuracy
np.mean(target == pred)
Out[15]:
0.7153587189249103
Parallélisation avec joblib (process)¶
In [16]:
from joblib import Parallel, delayed
#one_vs_rest parallélisée - technologie processus
def one_vs_rest_jb_process(Z, cible, max_iter=10, tol=1e-6):
nb_classes = np.max(cible)
#instanciation de la classe pour le traitement parallèle
#on peut restreindre le nombre de ressources à solliciter (n_jobs)
#technologie PROCESSUS
processes = Parallel(n_jobs=nb_classes,prefer="processes")
#le résultat est récupéré sous la forme d'une liste de vecteurs de coefficients
#les items sont dans l'ordre des soumissions, pas de réorganisation à faire
lst_coefs = processes(delayed(fit_reg)(Z,cible==(k+1),max_iter,tol) for k in range(nb_classes))
#return liste des coefficients
return lst_coefs
In [17]:
#entraînement
coefs_bis_ = one_vs_rest_jb_process(Z,target)
In [18]:
#vérification
np.mean(target == pred_reg(Z,np.array(coefs_bis_)))
Out[18]:
0.7153587189249103
Parallélisation avec joblib (thread)¶
In [19]:
#one_vs_rest parallélisée, technologie thread
def one_vs_rest_jb_thread(Z, cible, max_iter=10, tol=1e-6):
nb_classes = np.max(cible)
#instanciation de la classe pour le traitement parallèle
#technologie THREADS
threads = Parallel(n_jobs=nb_classes,prefer="threads")
#le résultat est récupéré sous la forme d'une liste cf. le param. "return_as"
#liste de vecteurs de coefficients ici
#les vecteurs sont dans l'ordre des soumissions et non pas d'achèvement
#par conséquent pas de réorganisation à faire
lst_coefs = threads(delayed(fit_reg)(Z,cible==(k+1),max_iter,tol) for k in range(nb_classes))
#return liste des coefficients
return lst_coefs
In [20]:
#entraînement
coefs_ter_ = one_vs_rest_jb_thread(Z,target)
In [21]:
#vérification
np.mean(target == pred_reg(Z,np.array(coefs_ter_)))
Out[21]:
0.7153587189249103
Durées d'exécution¶
In [22]:
# durée d'excution - SANS parallélisation
%timeit -r 10 -n 1 one_vs_rest(Z, target)
11.9 s ± 800 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)
In [23]:
# durée d'excution - AVEC parallélisation / PROCESSUS
%timeit -r 10 -n 1 one_vs_rest_jb_process(Z, target)
8.4 s ± 130 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)
In [24]:
# durée d'excution - AVEC parallélisation / THREAD
%timeit -r 10 -n 1 one_vs_rest_jb_thread(Z, target)
8.63 s ± 220 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)