Fondamentaux du Machine Learning Partie 2 : Régression Linéaire MultivariéeDraft

20 septembre 2023

Bienvenue dans la Partie 2 / N de ma série sur les Fondamentaux du Machine Learning en C ! Si ce n’est déjà fait, parcourez la partie 1 — je vais supposer que vous êtes familier avec les concepts présentés là-bas dans cet article.

Au-delà d’une seule variable

Notre premier modèle d’apprentissage automatique, le modèle de régression linéaire univariée, était techniquement de l’apprentissage automatique, mais pour être honnêtes, il n’est pas très impressionnant. Alors, quel était l’intérêt de tous ces efforts ? Jeter les bases pour des modèles plus complexes, comme le modèle linéaire multivariable que nous allons examiner aujourd’hui.

Le modèle linéaire multivarié

Le modèle linéaire multivarié (MLM) est essentiellement le même que son analogue univarié, à ceci près qu’il nous permet d’avoir plusieurs variables d’entrée qui influencent la sortie. Dans la dernière partie, nous avons essayé de prédire les prix des logements en nous basant uniquement sur la superficie. Nous savons tous qu’une maison de 5000 $pi^{2}$ en Arkansas est plusieurs fois moins chère qu’une maison similaire à San José, mais notre modèle donnerait le même prix pour les deux.

C’est là que nous réalisons que notre modèle doit être capable de combiner des variables de différentes natures, chacune pouvant avoir un effet différent sur notre prédiction.

Énoncé du problème

Commençons par définir notre problème :

On nous donne des données pour les variables d’entrée suivantes et leurs prix de maison correspondants $y_{n}$ (en milliers de dollars) :

$x_{1}$ : Surface habitable (en $pi^{2}$ )

$x_{2}$ : Largeur du terrain (en pi)

$x_{3}$ : Profondeur du terrain (en pi)

$x_{4}$ : Distance de la côte (en miles)

Prédire le prix d’une maison à partir d’un nouvel ensemble de données.

Nous savons que chacun de ces facteurs aura un effet sur le prix. Certains (comme la surface habitable et la distance de la côte) auront un effet très important, tandis que d’autres en auront un plus faible.

Alors, comment pouvons-nous exprimer ces relations sous forme de modèle ? Essayons d’en construire un en examinant quelques graphiques et en faisant des estimations approximatives.

Graphique : Surface habitable vs. Prix

On peut observer une association positive entre la surface habitable et le prix. La droite de tendance est modélisée par $0.34 x_{1} + 550$ .

Graphique : Largeur du terrain vs. Prix

On observe également ici une légère association positive. La droite de tendance est $3.3 x_{2} + 1771$ .

Graphique : Profondeur du terrain vs Prix

Cette relation positive est bien plus forte que pour la largeur du terrain, ce qui suggère qu’une forme de terrain plus profonde est davantage valorisée par les acheteurs. La droite de tendance est $5.84 x_{3} + 1459$ .

Graphique : Distance par rapport au rivage vs Prix

Comme on pouvait s’y attendre, plus on s’éloigne du rivage, plus le prix est bas. On constate également que cet effet diminue à mesure que l’on s’éloigne davantage. La droite de tendance est $- 8.07 x_{4} + 2758$ .

Une première estimation du modèle

Les lignes de tendance nous indiquent la relation approximative entre chaque caractéristique et le prix du logement. Nous pouvons les combiner pour obtenir un modèle approximatif des données. Essayons d’additionner toutes les lignes de tendance et voyons si nous obtenons un modèle raisonnable.

\overset{y}{^} = 0.34 x_{1} + 550 + 3.3 x_{2} + 1771 + 5.84 x_{3} + 1459 - 8.07 x_{4} + 2758 = 0.34 x_{1} + 3.3 x_{2} + 5.84 x_{3} - 8.07 x_4 + 6538

Intuitivement, qu’est-ce que cela signifie ? Nous avons une valeur de base (biais) de 6,5 M$, qui représente la valeur théorique d’une maison avec 0 m², sans terrain, au bord de l’océan. Attendez… cela ne semble pas correct. Nous devrions probablement prendre la moyenne des ordonnées à l’origine car chaque graphique se superpose aux autres. Nouveau modèle :

\overset{y}{^} = 0.34 x_{1} + 3.3 x_{2} + 5.84 x_{3} - 8.07 x_{4} + (550 + 1771 + 1459 + 2758) /4 = 0.34 x_{1} + 3.3 x_{2} + 5.84 x_{3} - 8.07 x_{4} + 1634.5

Ok… Maintenant c’est 1,6 M$. Mais c’est peut-être raisonnable puisqu’elle est juste au bord de l’eau. Le reste du modèle semble logique. Il suggère que

Pour chaque augmentation de 1 unité de	le prix du logement change de
Surface habitable	340 $
Largeur du terrain	3300 $
Profondeur du terrain	5840 $
Distance du rivage	-8070 $

Bien que raisonnables, ce ne sont que des hypothèses basées sur des calculs approximatifs. Pour véritablement résoudre notre problème, nous devons trouver le MLM optimal qui prédit les prix des logements à partir de ces 4 caractéristiques. En termes mathématiques, nous devons trouver les poids $w_{1}, w_{2}, w_{3}, w_{4}$ et le biais $b$ tels que

\overset{y}{^} = w_{1} x_{1} + w_{2} x_{2} + w_{3} x_{3} + w_{4} x_{4} + b

prédise les prix des logements avec une erreur minimale.

MLM, en code

Comme dans la partie 1, définissons notre modèle. Puisque nous avons plus d’un poids, nous devons utiliser un tableau de w.

struct mlinear_model {
	int num_weights;
	double *w;
	double b;
};

Maintenant, vous avez peut-être pu anticiper de quelques étapes et réaliser que si nous représentons w comme un vecteur, et x comme un vecteur, nous pouvons représenter la sortie du modèle comme le produit scalaire de $w$ et $x$ , additionné avec $b$

$$ \begin{align*} \hat{y} &= w_1 x_1 + \dots + w_nx_n + b \\ &= \vec{x} \cdot \vec{w} + b \\ \end{align*}


Ahh, beaucoup plus propre, non ? Enfin, pas pour longtemps...

Je vais faire un autre saut ici, et représenter notre vecteur de longueur $n$ comme une matrice $(n \times 1)$.

$$
w = \begin{bmatrix} w_1 \\\\ w_2 \\\\ w_3 \\\\ \vdots \\\\ w_n \end{bmatrix}

Vous verrez bientôt comment cela est utile. Définissons notre struct matrix.

// matrix.h

// Juste pour pouvoir remplacer `double` par n'importe quel flottant plus tard
typedef double mfloat;

typedef struct {
	// row major
	mfloat *buf;
	int rows;
	int cols;
} matrix;

Maintenant, le modèle ressemble à

// main.c
struct mlinear_model {
	matrix w;
	mfloat b;
};

Notez que num_weights est maintenant stocké dans la matrice.

Puisque nous sommes passés aux matrices, nous devons utiliser le produit matriciel au lieu du produit scalaire vectoriel. J’ai ajouté quelques notes qui expliquent le code, au cas où vous auriez besoin de rafraîchir votre algèbre linéaire.

// matrix.h

// Obtenir un élément dans la matrice
mfloat
matrix_get(matrix m, int row, int col)
{
    return m.buf[row * m.cols + col];
}
// Définir un élément dans la matrice
void
matrix_set(matrix m, int row, int col, mfloat val)
{
    m.buf[row * m.cols + col] = val;
}

// out = m1 dot m2
void
matrix_dot(matrix out, const matrix m1, const matrix m2)
{
	// Sur la i-ème ligne de la première matrice
    for (int row = 0; row < m1.rows; row++) {
		// Sur la j-ème colonne de la deuxième matrice
        for (int col = 0; col < m2.cols; col++) {
			// Exécuter le produit scalaire vectoriel et mettre le résultat dans out[i][j]
            mfloat sum = 0.0;
			// m1.cols == m2.rows donc k itère sur tout
            for (int k = 0; k < m1.cols; k++) {
                mfloat x1 = matrix_get(m1, row, k);
                mfloat x2 = matrix_get(m2, k, col);
                sum += x1 * x2;
            }
            matrix_set(out, row, col, sum);
        }
    }
}

Premièrement, notons les propriétés suivantes concernant le produit matriciel entre $X$ et $W$

Il ne peut être calculé que si et seulement si X.cols == W.rows
La matrice résultante a les dimensions (X.rows, W.cols)

Dans ce cas,

(1 \times n) X \cdot (n \times 1) \cdot W \to \to (1 \times 1) r \overset{e}{ˊ} sultat

Ainsi, $X$ doit être un vecteur ligne et $W$ doit être un vecteur colonne pour que nous puissions reproduire le comportement du produit scalaire vectoriel.

$$ \begin{align*} X \cdot W &= \begin{bmatrix} x_1 & x_2 & x_3 & \cdots & x_n \end{bmatrix} \cdot \begin{bmatrix} w_1 \\ w_2 \\ w_3 \\ \vdots \\ w_n \end{bmatrix} \\ &= x_1 w_1 + x_2 w_2 + \dots + x_nw_n \end{align*}


> Vous pouvez penser au produit matriciel comme prenant chaque colonne de la deuxième matrice, la retournant sur le côté, et faisant son produit scalaire avec la ligne de la première matrice.
>
> Jetez un œil au code et vérifiez que la multiplication ci-dessus renverra le résultat annoncé.

Maintenant, nous avons les outils pour coder la prédiction.

```c
// mlr.c

// x est un vecteur ligne, ou une matrice (1 x n)
mfloat predict(struct mlinear_model model, const matrix x) {
	// (1 x n) . (n x 1) => (1 x 1) matrice, qui est un nombre
	mfloat result[1][1] = {0.0};
	matrix tmp = {.buf = result, .rows = 1, .cols = 1};
	// Définir tmp comme le résultat
	matrix_dot(tmp, x, model.w);
	return tmp.buf[0] + model.b;
}

## Optimisation du modèle

On nous donne 2 tableaux : les données d'entrée, formatées comme une matrice $(m \times 4)$.

$$
X = \begin{bmatrix} x_1^1 & x_2^1 & x_3^1 & x_4^1 \\\\ x_1^2 & x_2^2 & x_3^2 & x_4^2 \\\\ \vdots & \vdots & \vdots & \vdots \\\\ x_1^m & x_2^m & x_3^m & x_4^m \end{bmatrix}
$$

Chaque ligne représente un échantillon, et les valeurs dans chaque colonne sont les caractéristiques décrites précédemment.

Et les prix des maisons correspondants $Y$, en milliers de dollars :

$$
Y = \begin{bmatrix} y^1 \\\\ y^2 \\\\ y^3 \\\\ y^4 \\\\ \vdots \\\\ y^m \end{bmatrix}
$$
Nous avons 100 échantillons de données, donc $m = 100$. Je vais les stocker dans `data.h`

```c
// data.h

#define NUM_FEATURES 4
#define NUM_SAMPLES 100

static mfloat X_data[NUM_SAMPLES][NUM_FEATURES] = { /* données omises*/ };

static mfloat Y_data[NUM_SAMPLES] = { /* données omises */ };

static matrix X = {.buf = (mfloat *)X_data,
                   .rows = NUM_SAMPLES,
                   .cols = NUM_FEATURES};

static matrix Y = {
    .buf = (mfloat *)Y_data, .rows = NUM_SAMPLES, .cols = 1};

Maintenant, c’est l’heure de l’optimisation !

Cela signifie que nous devons trouver le modèle avec les paramètres $W$ et $b$ tel que l’erreur sur tous les échantillons de données soit minimisée. Mais quelle est notre erreur ? Nous pouvons en fait utiliser la même définition que dans la partie 1, puisque nous comparons toujours deux nombres $\overset{y}{^}^{(i)}$ et $y^{(i)}$ .

J_{Wb} (X) = \frac{1}{m} i = 0 \sum m - 1 (\overset{y}{^}^{(i)} - y^{(i)})^{2}

C’est la moyenne de la somme des différences au carré, ou la différence quadratique moyenne entre les valeurs $y$ attendues et réelles. Plus cette valeur est proche de $0$ , meilleur est notre modèle. Réécrivons $\overset{y}{^}$ en fonction de $W, b,$ et $X$

J_{Wb} (X) = \frac{1}{m} i = 0 \sum m - 1 (X^{(i)} \cdot W + b - y^{(i)})^{2}

Ici, $X^{(i)}$ fait référence à la $i$ -ème ligne de $X$ . $X^{(i)} \cdot W$ est le produit scalaire du $i$ -ème échantillon avec les poids. Après avoir ajouté le biais $b$ , nous obtenons $\overset{y}{^}$ , ou la prédiction.

En code :

// Retourne une matrice d'une seule ligne qui représente la ligne i de m
matrix matrix_get_row(matrix m, int i) {
  return (matrix){
      .buf = m.buf + i * m.cols, .rows = 1, .cols = m.cols};
}

mfloat cost(struct mlinear_model *model, matrix X,
                    matrix Y) {
  mfloat cost = 0.0;
  for (int i = 0; i < X.rows; i++) {
    matrix x_i = matrix_get_row(X, i);
    mfloat f_wb = predict(model, x_i);
    mfloat diff = matrix_get(Y, 0, i) - f_wb;
    cost += diff * diff;
  }
  return cost / (2.0 * X.rows);
}

#### Quelle était la qualité de l'estimation ?

Revenons à notre estimation et voyons comment elle se comporte.

```c
int main() {
  matrix W = matrix_new(X.cols, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};
  printf("Coût du modèle nul : %f\n", cost(&model, X, Y));
  matrix_set(W, 0, 0, 0.34);
  matrix_set(W, 1, 0, 3.3);
  matrix_set(W, 2, 0, 5.84);
  matrix_set(W, 3, 0, -8.07);
  model.b = 1634.5;
  printf("Coût du modèle estimé : %f\n",
         cost(&model, X, Y));
}

Sortie :

Coût du modèle nul :        3340683.781483
Coût du modèle estimé : 2641267.466911

Qui l’aurait cru ! Il fait légèrement mieux qu’une simple supposition. Retour à la descente de gradient.

Apportons le calcul différentiel. Nous devons dériver $J$ par rapport à $W$ et $b$ pour savoir dans quelle direction, et avec quelle magnitude, nous devons ajuster les poids et les biais pour réduire l’erreur.

\frac{\partial J}{\partial b} (X) = \frac{2}{m} \sum_i = 0^{m - 1} (X^{(i)} \cdot W + b - y^{(i)}) = \frac{2}{m} \sum_i = 0^{m - 1} (\overset{y}{^}^{(i)} - y^{(i)})

Assez simple, non ? Pas tant pour $W$ . Puisque $W$ est une matrice, chaque poids dans $W$ a un effet distinct sur la sortie. Cela signifie que nous devons calculer une dérivée distincte pour chaque poids $W^{(j)}$ .

\frac{\partial J}{\partial W ^{(j)}} (X) = \frac{2}{m} i = 0 \sum m - 1 (X^{(i)} \cdot W + b - y^{(i)}) X_{j}^{(i)}

Ce terme $X_{j}^{(i)}$ provient de la dérivation de $X^{(i)} \cdot W$ . Pour un poids donné $W^{(j)}$ , et une ligne $i$ , il y aurait un terme de la forme $X_{j}^{(i)} W^{(j)}$ . Puisque $X_{j}^{(i)}$ est un coefficient constant dans la dérivée partielle par rapport à $W^{(j)}$ , il est extrait en raison de la règle de la chaîne.

Puisque nous pensons déjà en termes de matrices, écrivons la dérivée de $J$ par rapport à la matrice $W$ entière.

\frac{\partial J}{\partial W} = \frac{\partial J}{\partial W ^{(1)}} \frac{\partial J}{\partial W ^{(2)}} ⋮ \frac{\partial J}{\partial W ^{(n)}} = \frac{2}{m} \sum_i = 0^{m - 1} (X^{(i)} \cdot W + b - y^{(i)}) X^{(i)}_1 \frac{2}{m} \sum_i = 0^{m - 1} (X^{(i)} \cdot W + b - y^{(i)}) X^{(i)}_2 ⋮ \frac{2}{m} \sum_i = 0^{m - 1} (X^{(i)} \cdot W + b - y^{(i)}) X^{(i)}_n = \frac{2}{m} i = 0 \sum m - 1 (X^{(i)} \cdot W + b - y^{(i)}) X^{(i)}_1 X^{(i)}_2 ⋮ X^{(i)}_n = \frac{2}{m} \sum_i = 0^{m - 1} (X^{(i)} \cdot W + b - y^{(i)}) X^{(i)} = \frac{2}{m} \sum_i = 0^{m - 1} (\overset{y}{^}^{(i)} - y^{(i)}) X^{(i)}

Cela se simplifie bien ! Vous commencez peut-être à voir pourquoi nous aimons les matrices en apprentissage automatique. C’est une belle façon de penser aux calculs par lots.

Maintenant, pour le code

struct grad {
	matrix dJ_dW;
	mfloat dJ_db;
};

/* Calcule le gradient et écrit le résultat dans out. */
void compute_gradient(struct grad *out,
                      struct mlinear_model *model, const matrix X,
                      const matrix Y) {
  int m = X.rows;  // nombre d'échantillons
  int n = X.cols;  // nombre de caractéristiques

  // Utilisation de tmp pour stocker chaque ligne de X
  matrix tmp = matrix_new(1, n);
  for (int i = 0; i < m; i++) {
    // tmp = X^(i)
    matrix curr_row = matrix_get_row(X, i);
    // y_hat = (X^(i) dot W) + b
    mfloat y_hat = predict(model, curr_row);
    // yi = y^(i)
    mfloat yi = matrix_get(Y, 0, i);
    // Le terme entre parenthèses
    mfloat err = y_hat - yi;

    /*
     * Pour dJ_dW, nous devons multiplier l'erreur
     * par la ligne courante, et l'ajouter à la somme cumulée
     */

    // tmp = X^(i) * (y_hat^(i) - y^(i))
    matrix_scalar_mul(tmp, curr_row, err);
    // dJ_dW += tmp
    matrix_ip_T_add(out->dJ_dW, tmp);

    // dJ_db += (y_hat^(i) - y^(i))
    out->dJ_db += err;
  }

  /*
   * Je vais remplacer 2/m par 1/m ici puisque le facteur 2
   * peut être intégré à alpha dans l'étape suivante.
   */

  // dJ/db = (dJ/db) / m
  out->dJ_db /= m;
  // dJ/dW = (dJ/dW) / m
  matrix_scalar_ip_mul(out->dJ_dW, 1.0 / m);
  matrix_del(tmp);
}

Note : ip dans les fonctions de matrice signifie que le résultat est assigné au premier argument. Sinon, le résultat est assigné à un tampon passé en premier argument.

Les nouvelles fonctions matrix_* sont relativement simples, donc je vais omettre le code ici pour raccourcir l’article. En fait, ces fonctions ont été générées par des #define macros, donc le dépôt ne contient même pas le code source complet. Surveillez une Partie 2.5 qui les abordera !

Je vous encourage à parcourir le code et à vérifier qu’il correspond aux mathématiques.

Maintenant que nous avons le gradient, implémenter la descente de gradient est simple. Un rappel de l’algorithme de la partie 1 :

Répéter jusqu’à convergence :

Définir les valeurs initiales des poids et du biais : $W := [0 \dots 0]^{T}, b := 0$

Déplacer les variables dans la direction opposée aux dérivées, avec un pas $α$ : $W := W - α \frac{\partial J}{\partial W}, b := b - α \frac{\partial J}{\partial b}$

Aller à l’étape 2.

En code :

void gradient_descent(struct mlinear_model *model, const matrix X,
                      const matrix Y, const int num_iterations,
                      const mfloat alpha) {
  // tampon réutilisable pour le gradient
  int n = X.cols, m = X.rows;
  matrix dJ_dW = matrix_new(n, 1);
  struct grad tmp_grad = {.dJ_dW = dJ_dW, .dJ_db = 0.0};

  for (int i = 0; i < num_iterations; i++) {
    // Journalisation de la progression
    if (i % (num_iterations >> 4) == 0) {
      printf("\tCoût à l'itération %d : %f\n", i,
             compute_cost(model, X, Y));
    }
    // tmp_grad = gradient actuel pour le modèle
    compute_gradient(&tmp_grad, model, X, Y);
    // dJ/dW *= -alpha
    matrix_scalar_ip_mul(tmp_grad.dJ_dW, -alpha);
    // W += dJ/dW
    matrix_ip_add(model->W, tmp_grad.dJ_dW);
    // b += -alpha * dJ/db
    model->b += -alpha * tmp_grad.dJ_db;
  }
  matrix_del(dJ_dW);
}

Et c’est tout ! Maintenant, nous exécutons le programme sur les données pour voir comment notre coût s’améliore :

int main() {
  // hyperparamètres
  const int num_iterations = 1e7;
  const mfloat alpha = 1e-8;

  int n = X.cols, m = X.rows;
  matrix W = matrix_new(n, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};

  printf("Coût initial : %f\n", compute_cost(&model, X, Y));
  gradient_descent(&model, X, Y, num_iterations, alpha);
  printf("Coût final : %f\n", compute_cost(&model, X, Y));
  printf("Paramètres du modèle :\n");
  matrix_print(model.W);
  printf(" b=%f\n", model.b);
}

Sortie :

Coût initial : 3340683.781483
        Coût à l'itération 0 : 3340683.781483
        Coût à l'itération 625000 : 161369.409722
        Coût à l'itération 1250000 : 161332.630253
        Coût à l'itération 1875000 : 161315.326515
        Coût à l'itération 2500000 : 161298.041198
        Coût à l'itération 3125000 : 161280.768143
        Coût à l'itération 3750000 : 161263.507342
        Coût à l'itération 4375000 : 161246.258784
        Coût à l'itération 5000000 : 161229.022462
        Coût à l'itération 5625000 : 161211.798366
        Coût à l'itération 6250000 : 161194.586488
        Coût à l'itération 6875000 : 161177.386819
        Coût à l'itération 7500000 : 161160.199351
        Coût à l'itération 8125000 : 161143.024075
        Coût à l'itération 8750000 : 161125.860983
        Coût à l'itération 9375000 : 161108.710065
Coût final : 161091.571313
Paramètres du modèle :
W=[ 2.1185e-01 ]
[ 5.8627e+00 ]
[ 4.5904e+00 ]
[ -1.1702e+01 ]

 b=5.310395

Excellent ! Nous pouvons voir qu’à chaque itération notre coût diminue. Le modèle final ressemble à

\overset{y}{^} = .212 x_{1} + 5.86 x_{2} + 4.59 x_{3} - 11.7 x_{4} + 5.31

Optimiser l’optimiseur

Vous avez peut-être remarqué que même à la dix millionième itération de la descente de gradient, le coût continuait de diminuer. Pourquoi est-il si difficile de trouver les paramètres optimaux, surtout lorsque nous n’avons que 4 poids à optimiser ?

Une grande partie est due à la différence d’étalement des caractéristiques. Regardez ce diagramme en boîte

Nous pouvons voir que le poids pour la superficie doit se déplacer beaucoup plus que les poids pour les autres variables. Cependant, le pas d’apprentissage $α$ est le même lors de la mise à jour de chacun des poids. Cela signifie que nous devrons attendre longtemps après que les 3 autres poids de caractéristiques aient convergé pour que le poids de la superficie converge.

Une solution consiste à choisir des valeurs $α$ proportionnelles à l’étalement et à déplacer chaque poids par son $α$ correspondant. Mais cela signifierait stocker un autre tableau de valeurs alpha, ce qui est coûteux pour les grands modèles.

Une meilleure solution est de modifier nos données d’entrée pour qu’elles aient des étalements similaires. Nous pouvons le faire en les ajustant à une distribution normale. D’abord, nous calculons la moyenne de chaque caractéristique.

μ^{(i)} = \frac{1}{n} j = 1 \sum n X_{j}^{(i)}

Ensuite, nous calculons l’étalement, ou écart type, de chaque caractéristique.

σ^{(i)} = \frac{\sum _{j = 1}^{n} ( X _{j}^{(i)} - μ ^{(i)} ) ^{2}}{n}

Et enfin, nous la normalisons

Z^{(i)} = \frac{X ^{(i)} - μ ^{(i)}}{σ ^{(i)}}

Soustraire chaque colonne par sa moyenne centre les données afin qu’il y ait un étalement similaire de chaque côté de la distribution. Diviser par $σ$ rend l’étalement de chaque colonne approximativement le même. Écrivons cela en termes d’opérations matricielles

// Normaliser les données d'entrée `X`
void z_score_normalize(matrix X) {
  int n = X.cols, m = X.rows;
  // Tampon pour mu
  matrix mean = matrix_new(1, n);
  // Tampon pour sigma
  matrix stdev = matrix_new(1, n);

  // Calculer mu
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix_ip_add(mean, matrix_get_row(X, i));
  }
  matrix_scalar_ip_mul(mean, 1.0 / NUM_SAMPLES);

  // Calculer sigma
  matrix buf = matrix_new(1, n);
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix row = matrix_get_row(X, i);
    matrix_sub(buf, mean, row);
    matrix_ip_square(buf);
    matrix_ip_add(stdev, buf);
  }
  matrix_ip_sqrt(stdev);

  // Calculer Z
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix row = matrix_get_row(X, i);
    matrix_ip_sub(row, mean);
    matrix_ip_div(row, stdev);
  }
}

À quoi ressemblent nos données normalisées ?

Nous pouvons voir que les étalements sont maintenant tous à la même échelle. Comment cela améliore-t-il nos performances ?

int main() {
  int n = X.cols, m = X.rows;
  z_score_normalize(X);
  matrix W = matrix_new(n, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};
  printf("Coût initial : %f\n", compute_cost(&model, X, Y));
  const int num_iterations = 1e7;
  const mfloat alpha = 1e-8;
  gradient_descent(&model, X, Y, num_iterations, alpha);
  printf("Coût final : %f\n", compute_cost(&model, X, Y));

  printf("Paramètres du modèle :\nW=");
  matrix_print(model.W);
  printf(" b=%f\n", model.b);
}

Sortie :

Coût initial : 3340683.781483
        Coût à l'itération 0 : 3340683.781483
        Coût à l'itération 625000 : 3305547.933747
        Coût à l'itération 1250000 : 3270852.031317
        Coût à l'itération 1875000 : 3236590.554988
        Coût à l'itération 2500000 : 3202758.054240
        Coût à l'itération 3125000 : 3169349.146943
        Coût à l'itération 3750000 : 3136358.518493
        Coût à l'itération 4375000 : 3103780.920969
        Coût à l'itération 5000000 : 3071611.172296
        Coût à l'itération 5625000 : 3039844.155415
        Coût à l'itération 6250000 : 3008474.817473
        Coût à l'itération 6875000 : 2977498.169013
        Coût à l'itération 7500000 : 2946909.283181
        Coût à l'itération 8125000 : 2916703.294939
        Coût à l'itération 8750000 : 2886875.400289
        Coût à l'itération 9375000 : 2857420.855511
Coût final : 2828334.976402
Paramètres du modèle :
W=[ 8.3053e+00 ]
[ 2.4385e+00 ]
[ 6.0186e+00 ]
[ -2.2844e+00 ]

 b=227.136534

C’est encore plus lent !? Mais attendez—regardez ce qui se passe lorsque nous changeons $α = 1 0^{- 8}$ en $α = 1,$ ce qui aurait fait diverger notre ancien modèle.

Coût initial : 3340683.781483
        Coût à l'itération 0 : 3340683.781483
        Coût à l'itération 625000 : 136947.544865
        Coût à l'itération 1250000 : 136947.544865
        Coût à l'itération 1875000 : 136947.544865
        Coût à l'itération 2500000 : 136947.544865
        Coût à l'itération 3125000 : 136947.544865
        Coût à l'itération 3750000 : 136947.544865
        Coût à l'itération 4375000 : 136947.544865
        Coût à l'itération 5000000 : 136947.544865
        Coût à l'itération 5625000 : 136947.544865
        Coût à l'itération 6250000 : 136947.544865
        Coût à l'itération 6875000 : 136947.544865
        Coût à l'itération 7500000 : 136947.544865
        Coût à l'itération 8125000 : 136947.544865
        Coût à l'itération 8750000 : 136947.544865
        Coût à l'itération 9375000 : 136947.544865
Coût final : 136947.544865
Paramètres du modèle :
W=[ 8.1088e+03 ]
[ 8.1764e+02 ]
[ 6.6840e+02 ]
[ -3.6835e+03 ]

 b=2364.131545

Wow ! Dès le premier log, il a convergé vers la valeur minimale ! En fait, avec $α = 1$ , nous n’avons besoin que d’environ $3 \times 1 0^{4}$ itérations pour la convergence au lieu de $1 0^{7}$ ! Nous pouvons donc voir que l’utilisation de la normalisation z-score peut accélérer la descente de gradient de plusieurs ordres de grandeur, surtout lorsque les caractéristiques ont des étalements variables.

Si vous avez remarqué que les nouveaux $W$ et $b$ sont très éloignés des originaux, c’est parce que nous avons fondamentalement changé les données d’entrée. Nous modélisons maintenant basé sur l’amplitude de l’écart d’un échantillon de caractéristique au lieu de sa valeur absolue.

Conclusion

Voilà donc tout pour la régression linéaire multivariée. Ce sera de loin l’article le plus difficile de la série, alors n’ayez crainte pour la suite ! Tout le code est disponible ici. Je vous encourage à l’exécuter et à en modifier le comportement.

Si vous avez des questions ou des commentaires, n’hésitez pas à laisser un commentaire ou à m’envoyer un e-mail.

✦ Aucune IA n'a été utilisée dans la conception, la recherche, la rédaction ou l'édition de cet article.