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]:
#scikit-learn
import sklearn
sklearn.__version__
Out[3]:
'1.5.1'
In [4]:
#numba
import numba
numba.__version__
Out[4]:
'0.60.0'

Chargement et préparation des données¶

In [5]:
import sklearn.datasets as ds
covtype = ds.fetch_covtype()
In [6]:
#matrice des descripteurs
X = np.array(covtype.data)
X.shape
Out[6]:
(581012, 54)
In [7]:
#attribut classe
target = np.array(covtype.target)
target.shape
Out[7]:
(581012,)
In [8]:
#valeurs possibles
np.unique(target,return_counts=True)
Out[8]:
(array([1, 2, 3, 4, 5, 6, 7]),
 array([211840, 283301,  35754,   2747,   9493,  17367,  20510],
       dtype=int64))
In [9]:
#standardisation des X
from sklearn.preprocessing import StandardScaler
std = StandardScaler()
Z = std.fit_transform(X)

#dim.
Z.shape
Out[9]:
(581012, 54)

Scikit-Learn¶

In [10]:
from sklearn.linear_model import LogisticRegression
from sklearn.multiclass import OneVsRestClassifier
#instanciation - n_jobs = 1 -> 1 seul coeur sollicité, pas de parallélisation
ove_first = OneVsRestClassifier(LogisticRegression(random_state=0), n_jobs = 1)
#entraînement
ove_first.fit(Z,target)
#vérif. liste des classes
ove_first.classes_
Out[10]:
array([1, 2, 3, 4, 5, 6, 7])
In [11]:
#accuracy
np.mean(target == ove_first.predict(Z))
Out[11]:
0.7153724880036901
In [12]:
#avec parallélisation - n_jobs = -1 -> tous les coeurs sollicités
ove_bis = OneVsRestClassifier(LogisticRegression(random_state=0), n_jobs = -1)
#entraînement
ove_bis.fit(Z,target)
#accuracy
np.mean(target == ove_bis.predict(Z))
Out[12]:
0.7153724880036901
In [13]:
# durée d'exécution - SANS parallélisation
%timeit -r 7 -n 1 ove_first.fit(Z,target)
14.3 s ± 713 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [14]:
# durée d'exécution - AVEC parallélisation
%timeit -r 7 -n 1 ove_bis.fit(Z,target)
11.8 s ± 519 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

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 (Numba)¶

In [15]:
from numba import njit

@njit
def fit_reg_numba(Z, y, max_iter=5, 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 [16]:
#on part du principe que les classes sont identifiées par 1, 2, 3, ....
@njit
def one_vs_rest(Z, cible, max_iter=5, 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_numba(Z,y,max_iter,tol))
    #return liste des coefficients
    return lst_coefs
In [17]:
#entraînement sur nos données
coefs_ = one_vs_rest(Z,target)
print(len(coefs_))
7
In [18]:
#sous-forme matricielle
betas_ = np.array(coefs_)
betas_.shape
Out[18]:
(7, 55)

Fonction de prédiction¶

In [19]:
# 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 [20]:
# appliquer
pred = pred_reg(Z,betas_)
np.unique(pred,return_counts=True)
Out[20]:
(array([1, 2, 3, 4, 5, 6, 7], dtype=int64),
 array([205285, 305670,  50098,    675,    393,   3264,  15627],
       dtype=int64))
In [21]:
# accuracy
np.mean(target == pred)
Out[21]:
0.7143845566012406

Numba avec parallélisation¶

In [22]:
#prange pour la paralléisation des boucles
from numba import prange
#modification de la directive
@njit(parallel=True)
def one_vs_rest_parallel(Z, cible, max_iter=5, tol=1e-6):
    nb_classes = np.max(cible)
    #pour chaque modalité de la variable cible
    #créer une liste de vecteurs numpy
    lst_coefs = [np.empty(Z.shape[1]+1)] * nb_classes
    for k in prange(nb_classes):
        y = (cible==(k+1))
        #remplacer les vecteurs par ceux estimés
        lst_coefs[k] = fit_reg_numba(Z,y,max_iter,tol)
    #return liste des coefficients
    return lst_coefs
In [23]:
#entraînement
coefs_bis_ = one_vs_rest_parallel(Z,target)
In [24]:
#vérification
np.mean(target == pred_reg(Z,np.array(coefs_bis_)))
Out[24]:
0.7143845566012406
In [26]:
# durée d'excution - SANS parallélisation
%timeit -r 7 -n 1 one_vs_rest(Z, target)
17.4 s ± 142 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [27]:
# durée d'excution - AVEC parallélisation
%timeit -r 7 -n 1 one_vs_rest_parallel(Z, target)
9.3 s ± 694 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)