Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Le perceptron est la brique élémentaire des réseaux de neurones. Il calcule une somme pondérée de ses entrées, ajoute un biais, et passe le résultat dans une fonction d’activation. Dans ce carnet, nous :

  1. Construisons un perceptron à partir de zéro avec NumPy.
  2. L’entraînons avec la règle d’apprentissage classique du perceptron pour séparer des événements sismiques du bruit.
  3. Réutilisons le même perceptron comme régresseur linéaire entraîné par descente de gradient, et le comparons aux moindres carrés ordinaires (ordinary least squares, OLS).
  4. Étudions comment le taux d’apprentissage (learning rate) modifie l’entraînement.

🖥️ Diapositives du cours — Séance 21 (mer. 18 nov.)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

import mlgeo_synth

1. Un perceptron à partir de zéro

Un perceptron d’entrées x1,,xnx_1, \dots, x_n, de poids w1,,wnw_1, \dots, w_n et de poids de biais bb calcule

o=ϕ(jwjxj+b),o = \phi\left(\sum_j w_j x_j + b\right),

ϕ est la fonction d’activation. Nous implémentons trois activations : aucune (identité), sigmoïde et échelon.

class simplePerceptron:
    # Activation functions
    def __none(x):
        return x
    def __sigmoid(x):
        return 1 / (1 + np.exp(-x))
    def __step(x):
        out = np.zeros(x.shape)
        out[x >= 0] = 1
        return out
    possibleActivations = {
        'none': __none,
        'sigmoid': __sigmoid,
        'step': __step
    }
    # Initialize
    def __init__(self, w=None, b_w=None, activation=None):
        if w is None:
            self.w = 0
        else:
            self.w = w
        if b_w is None:
            self.b_w = 0
        else:
            self.b_w = b_w
        if (activation is None) or (activation not in self.possibleActivations):
            self.activation = 'none'
        else:
            self.activation = activation
    def setWs(self, w, b_w):
        self.w = w
        self.b_w = b_w
    def predict(self, x):
        # x should be an (m, n) array: each row is one observation
        weightedSum = np.zeros(x.shape[0])
        for i in range(len(self.w)):
            weightedSum += self.w[i] * x[:, i]
        summed = weightedSum + (1 * self.b_w)
        return self.possibleActivations[self.activation](summed)

1.1 Tester le perceptron

Avec un poids fixé à 1, un biais de 0 et aucune activation, le perceptron doit retourner son entrée inchangée.

ourPerceptron = simplePerceptron(w=[0], b_w=1, activation='none')
ourPerceptron.setWs([1], 0)
testInput = np.array([[100, -1, 0.1]]).T
ourPerceptron.predict(testInput)
array([100. , -1. , 0.1])

2. Classifier des détections sismiques

Les détecteurs sismiques résument de courtes fenêtres d’un sismogramme par quelques caractéristiques (features). La table mlgeo_synth.detector_features contient quatre caractéristiques par fenêtre : le rapport STA/LTA (sta_lta), le kurtosis de la forme d’onde, le centroïde spectral en Hz et la fréquence dominante en Hz. L’étiquette vaut 1 pour une fenêtre de séisme et 0 pour une fenêtre de bruit, avec des classes équilibrées.

Nous utilisons deux caractéristiques, sta_lta et spectral_centroid_hz. Les événements ont des arrivées impulsives (STA/LTA élevé) et plus d’énergie haute fréquence (centroïde spectral élevé) que le bruit de fond.

df = mlgeo_synth.detector_features(n=2000, event_fraction=0.5, seed=0)
print(df.head())

fig, ax = plt.subplots(figsize=(6, 4.5))
for lbl, name in [(0, 'noise'), (1, 'event')]:
    part = df[df['label'] == lbl]
    ax.scatter(part['sta_lta'], part['spectral_centroid_hz'], s=8, alpha=0.5, label=name)
ax.set_xlabel('STA/LTA')
ax.set_ylabel('Spectral centroid (Hz)')
ax.set_title('Detector features: all 2000 windows')
ax.legend();
     sta_lta   kurtosis  spectral_centroid_hz  dominant_freq_hz  label
0  10.419133  22.113184              9.556484          7.786286      1
1   3.314406   2.773164              1.716694          1.009758      0
2   2.336543   4.639630              1.271397          0.251278      0
3   3.976867  13.591792              6.692700          4.682170      1
4   1.151811   3.131239              4.474578          1.426544      0
<Figure size 600x450 with 1 Axes>

2.1 Pourquoi une tranche facile à deux caractéristiques ?

La règle d’apprentissage classique du perceptron ne converge que lorsque les données sont linéairement séparables : il doit exister une droite qui classe correctement chaque échantillon d’entraînement. Si les classes se recouvrent, la règle continue de mettre à jour les poids indéfiniment. Le nuage de points ci-dessus montre une petite bande de recouvrement entre les deux nuages : le jeu de données complet n’est donc pas linéairement séparable.

Nous sélectionnons donc un sous-ensemble quasi linéairement séparable : nous calculons un score combiné simple, sta_lta + spectral_centroid_hz, et écartons les échantillons de la bande de recouvrement. Nous nous en tenons aussi à exactement deux caractéristiques, pour pouvoir tracer la frontière de décision comme une droite dans le plan. Les données réelles de détecteurs sont plus désordonnées ; les carnets suivants traitent le recouvrement avec des réseaux multicouches et des pertes fondées sur le gradient.

score = df['sta_lta'] + df['spectral_centroid_hz']
keep = ((df['label'] == 0) & (score < 5)) | ((df['label'] == 1) & (score > 9))
subset = df[keep]

# Balance the classes: 300 windows of each
rng = np.random.default_rng(42)
noise_idx = rng.choice(subset.index[subset['label'] == 0], 300, replace=False)
event_idx = rng.choice(subset.index[subset['label'] == 1], 300, replace=False)
subset = subset.loc[np.concatenate([noise_idx, event_idx])]

simpleInputs = subset[['sta_lta', 'spectral_centroid_hz']].to_numpy()
simpleOutputs = subset['label'].to_numpy()

# Shuffle the samples
p = rng.permutation(len(simpleInputs))
simpleInputs = simpleInputs[p, :]
simpleOutputs = simpleOutputs[p]

fig, ax = plt.subplots(figsize=(6, 4.5))
for lbl, name in [(0, 'noise'), (1, 'event')]:
    m = simpleOutputs == lbl
    ax.scatter(simpleInputs[m, 0], simpleInputs[m, 1], s=8, alpha=0.6, label=name)
ax.set_xlabel('STA/LTA')
ax.set_ylabel('Spectral centroid (Hz)')
ax.set_title('Balanced, near-linearly-separable subset (600 windows)')
ax.legend();
<Figure size 600x450 with 1 Axes>

2.2 La règle d’apprentissage du perceptron

Nous partons de poids tous nuls et utilisons l’activation en échelon, de sorte que le perceptron sort 0 ou 1. Nous balayons ensuite le jeu de données. Pour chaque échantillon, nous :

  1. Faisons une prédiction avec les poids courants.

  2. Mettons à jour chaque poids selon la règle

    wj=wj+Δwjw_j = w_j + \Delta w_j, où Δwj=η(tioi)xji\Delta w_j = \eta\,(t^i - o^i)\,x^i_j

    Ici η est le taux d’apprentissage, tit^i et oio^i sont la cible et la sortie pour l’échantillon ii, et xjix^i_j est l’entrée jj de cet échantillon. Le poids de biais reçoit la même mise à jour, avec une entrée égale à 1. Notez que les échantillons correctement classés (ti=oit^i = o^i) laissent les poids inchangés.

  3. Répétons des balayages complets (des époques) sur les données jusqu’à ce qu’une époque ne produise aucune erreur.

Un détail : avec des poids initiaux nuls, le taux d’apprentissage ne fait que mettre à l’échelle le vecteur de poids final. L’activation en échelon ignore cette échelle : la suite des prédictions, et le nombre d’époques avant convergence, sont donc les mêmes pour tout η>0\eta > 0.

percept = simplePerceptron(w=np.zeros(2), b_w=0.0, activation='step')

learningRate = 0.01
maxEpochs = 50
mistakesPerEpoch = []

for epoch in range(maxEpochs):
    mistakes = 0
    for i in range(len(simpleInputs)):
        thisInput = simpleInputs[i:i + 1, :]
        thisTarget = simpleOutputs[i]
        out = percept.predict(thisInput)[0]
        error = thisTarget - out
        if error != 0:
            mistakes += 1
            percept.setWs(percept.w + learningRate * error * thisInput[0],
                          percept.b_w + learningRate * error)
    mistakesPerEpoch.append(mistakes)
    if mistakes == 0:
        break

accuracy = (percept.predict(simpleInputs) == simpleOutputs).mean()
print(f'Converged after {len(mistakesPerEpoch)} epochs; mistakes per epoch: {mistakesPerEpoch}')
print(f'Weights: {percept.w}, bias: {percept.b_w:.3f}, training accuracy: {accuracy:.3f}')
Converged after 2 epochs; mistakes per epoch: [42, 0]
Weights: [0.03415388 0.03796047], bias: -0.220, training accuracy: 1.000

La règle converge en quelques époques et classe correctement chaque échantillon d’entraînement. Comme nous utilisons deux caractéristiques, la frontière de décision est la droite w1x1+w2x2+b=0w_1 x_1 + w_2 x_2 + b = 0, que nous pouvons tracer directement.

fig, ax = plt.subplots(figsize=(6, 4.5))
for lbl, name in [(0, 'noise'), (1, 'event')]:
    m = simpleOutputs == lbl
    ax.scatter(simpleInputs[m, 0], simpleInputs[m, 1], s=8, alpha=0.6, label=name)

# Decision boundary: w1*x1 + w2*x2 + b = 0  ->  x2 = -(w1*x1 + b) / w2
x1 = np.linspace(simpleInputs[:, 0].min(), simpleInputs[:, 0].max(), 100)
x2 = -(percept.w[0] * x1 + percept.b_w) / percept.w[1]
ax.plot(x1, x2, 'k--', label='decision boundary')

ax.set_xlabel('STA/LTA')
ax.set_ylabel('Spectral centroid (Hz)')
ax.set_ylim(-0.5, simpleInputs[:, 1].max() + 0.5)
ax.set_title('Perceptron decision boundary')
ax.legend();
<Figure size 600x450 with 1 Axes>

Nous venons d’implémenter un algorithme en ligne (online) : les poids sont mis à jour après chaque échantillon, et non après avoir vu tout le jeu de données.

Exercice. Remplacez spectral_centroid_hz par kurtosis dans la cellule de sélection du sous-ensemble (gardez le même filtre par score sur les caractéristiques d’origine). La règle d’apprentissage converge-t-elle encore en moins de 50 époques ? Pourquoi, ou pourquoi pas ?

3. Ajuster une droite par descente de gradient

Le même perceptron, sans activation, est un modèle linéaire o=wx+bo = w x + b. Au lieu de la règle du perceptron, nous pouvons l’entraîner par descente de gradient sur une fonction de perte quadratique. Aucun terme de régularisation n’étant ajouté ici, cette perte moyennée sur les exemples est aussi la fonction objectif que minimise la descente de gradient. Nous générons un jeu de données linéaire bruité à ajuster.

num = 100
xRange = [0, 5]
yRange = [0, 2.5]

rng = np.random.default_rng(1)

x = np.linspace(xRange[0], xRange[1], num=num)
y = np.linspace(yRange[0], yRange[1], num=num)

# Add uniform noise in [-0.5, 0.5]
noise_factor = 0.5
x = x + (rng.random(num) * 2 - 1) * noise_factor
y = y + (rng.random(num) * 2 - 1) * noise_factor

# Shuffle x and y together
p = rng.permutation(num)
x = x[p]
y = y[p]
fig, ax = plt.subplots(figsize=(6, 4))
ax.scatter(x, y, color='k', s=12)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')
ax.set_axisbelow(True)
ax.grid(color='gray', linestyle='dashed')
<Figure size 600x400 with 1 Axes>

3.1 Fonction de perte

Nous utilisons l’erreur quadratique moyenne (mean squared error, MSE) comme perte, moyennée sur les exemples (mseCost ci-dessous).

def mseCost(prediction, target):
    mse = (np.square(prediction - target)).mean()
    return mse

# Sanity check: the cost of a perfect prediction is zero
mseCost(y, y)
np.float64(0.0)

3.2 Descente de gradient

La descente de gradient répète trois étapes : prédire, mesurer la perte, puis déplacer les poids d’un petit pas dans le sens opposé au gradient de la perte. Pour la MSE et un modèle linéaire, la mise à jour des poids est Δw=ηXT(to)\Delta w = \eta\, X^T (t - o) et celle du biais Δb=ηi(tioi)\Delta b = \eta \sum_i (t^i - o^i). Nous nous arrêtons quand la variation de la perte passe sous une tolérance, quand la perte cesse d’être finie (divergence), ou quand la limite d’itérations est atteinte.

def gradientDescent(perceptron, costFunction, trainInput, trainTarget,
                    learningRate, numIterations, stoppingCriterion):
    weights = []
    biases = []
    costs = []
    previousCost = None
    for i in range(numIterations):
        # Run the prediction
        prediction = perceptron.predict(trainInput)
        # Determine the cost
        thisCost = costFunction(prediction, trainTarget)
        # Stop if the cost diverged
        if not np.isfinite(thisCost):
            break
        # Stop if the change in cost is below the stopping criterion
        if previousCost and np.absolute(previousCost - thisCost) <= stoppingCriterion:
            break
        previousCost = thisCost
        # Record this weight, bias, and cost
        weights.append(perceptron.w)
        biases.append(perceptron.b_w)
        costs.append(thisCost)
        # Errors
        er = np.subtract(trainTarget, prediction)
        # Weight and bias updates
        weightUpdate = learningRate * np.dot(trainInput.T, er)
        biasWeightUpdate = learningRate * np.sum(er)
        perceptron.setWs(perceptron.w + weightUpdate, perceptron.b_w + biasWeightUpdate)
    return {'weights': weights, 'biases': biases, 'costs': costs}
numIterations = 10000
# Initialize the weight and bias to zero
ourPerceptron.setWs([0], 0)
# Train on the first 50 points, keep the rest for testing
output = gradientDescent(ourPerceptron, mseCost, np.reshape(x[0:50], (-1, 1)),
                         y[0:50], 0.001, numIterations, 1e-6)
# A function to plot outputs from gradient descent
def plotOutputs(gdOutput):
    iterationsRan = len(gdOutput['costs'])
    toPlot = ['weights', 'biases', 'costs']
    fig, axes = plt.subplots(1, 3, figsize=(12, 4))
    for i, key in enumerate(toPlot):
        axes[i].plot(range(iterationsRan), gdOutput[key], color='r')
        axes[i].set_xlabel('Iteration #')
        axes[i].set_ylabel('Value')
        axes[i].set_title(key.capitalize())
        axes[i].set_axisbelow(True)
        axes[i].xaxis.grid(color='gray', linestyle='dashed')
    fig.tight_layout()
    return fig

plotOutputs(output);
<Figure size 1200x400 with 3 Axes>
# Try the perceptron on data it has not seen before
testPredictions = ourPerceptron.predict(np.reshape(x[50:100], (-1, 1)))
testCost = mseCost(testPredictions, np.reshape(y[50:100], (-1, 1)))
'The mean squared error of our perceptron on the test half is: %f' % testCost
'The mean squared error of our perceptron on the test half is: 1.122684'

3.3 Comparaison avec les moindres carrés ordinaires

Les moindres carrés ordinaires (OLS) résolvent le même problème sous forme analytique. Un perceptron entraîné par descente de gradient, s’il fonctionne, doit aboutir à presque la même droite.

def compareOutputs(x, y, perceptron, xRange):
    fig, ax = plt.subplots(figsize=(12, 4))
    ax.scatter(x, y, color='k', s=12)
    ax.set_aspect('equal')
    forLine = np.linspace(xRange[0], xRange[1])
    ax.plot(forLine, perceptron.predict(np.reshape(forLine, (-1, 1))),
            color='r', linestyle='dotted', linewidth=1.5)
    # How well do we do relative to OLS?
    OLSoutput = LinearRegression().fit(x[0:50].reshape(-1, 1), y[0:50])
    ax.plot(forLine, OLSoutput.predict(forLine.reshape(-1, 1)),
            color='b', linestyle='dashed', linewidth=1.5)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_axisbelow(True)
    ax.grid(color='gray', linestyle='dashed')
    ax.legend(['Data', 'Perceptron fit', 'OLS fit'])
    return fig

compareOutputs(x, y, ourPerceptron, xRange);
<Figure size 1200x400 with 1 Axes>

4. Comment le taux d’apprentissage modifie l’entraînement

Le taux d’apprentissage η contrôle la taille du pas de la descente de gradient. Trop petit, l’entraînement se traîne ; trop grand, la perte oscille ou diverge. La grille ci-dessous entraîne le même perceptron, depuis le même départ à zéro, avec quatre taux d’apprentissage, et trace la courbe de perte de chacun (échelle logarithmique sur l’axe de la perte).

learningRates = [1e-5, 1e-4, 1e-3, 1e-2]

fig, axes = plt.subplots(2, 2, figsize=(10, 7))
with np.errstate(over='ignore', invalid='ignore'):
    for ax, lr in zip(axes.ravel(), learningRates):
        gdPerceptron = simplePerceptron(w=[0], b_w=0, activation='none')
        gdOut = gradientDescent(gdPerceptron, mseCost, np.reshape(x[0:50], (-1, 1)),
                                y[0:50], lr, 2000, 1e-8)
        costs = np.array(gdOut['costs'])
        # Truncate the curve once the cost exceeds 1e12 (divergence)
        plotCosts = costs[costs < 1e12]
        ax.semilogy(plotCosts, color='r', label=f'final cost = {costs[-1]:.3g}')
        ax.set_xlabel('Iteration #')
        ax.set_ylabel('MSE cost')
        ax.set_title(f'learning rate = {lr:g}')
        ax.legend()
        ax.set_axisbelow(True)
        ax.grid(color='gray', linestyle='dashed')
fig.tight_layout()
<Figure size 1000x700 with 4 Axes>

Lisez les quatre panneaux : à η=105\eta = 10^{-5}, la perte baisse encore après 2 000 itérations (trop lent). À η=104\eta = 10^{-4}, elle converge mais demande plus d’un millier d’itérations. À η=103\eta = 10^{-3}, elle converge en quelques centaines d’itérations. À η=102\eta = 10^{-2}, la perte croît sans borne : les pas dépassent le minimum et l’entraînement diverge (la courbe s’arrête là où la perte devient infinie).

4.1 Explorateur interactif

Le widget ci-dessous vous permet de faire varier le taux d’apprentissage, le nombre d’itérations, le critère d’arrêt et les poids initiaux, puis de réentraîner et de retracer. Les widgets ne fonctionnent pas dans la version statique du livre : téléchargez ce notebook et exécutez-le dans Jupyter pour l’utiliser. La grille statique ci-dessus montre la même leçon.

import ipywidgets as widgets
from IPython.display import display, clear_output
%matplotlib inline
# Create inputs
learningRateSlider = widgets.FloatLogSlider(
    value=.001,
    min=-6,
    max=1.0,
    step=1,
    description='Learning Rate:'
)

iterationsSlider = widgets.Dropdown(
    options=[1, 10, 100, 1000, 10000, 100000],
    value=100,
    description='Number of Iterations:'
)

stoppingCriterionSlider = widgets.FloatLogSlider(
    value=1e-4,
    min=-8,
    max=1,
    step=1,
    description='Stopping Criterion:'
)

startingWeight = widgets.BoundedFloatText(
    value=0,
    min=-50,
    max=50,
    description='Weight'
)

startingBias = widgets.BoundedFloatText(
    value=0,
    min=-50,
    max=5000,
    description='Bias'
)

# Create the "Update" button
updateBtn = widgets.Button(description="Update")

outputWidget = widgets.Output()

# Define the function to run when the "Update" button is clicked
def updateClick(_):
    with outputWidget:
        clear_output(wait=True)
        learningRate = learningRateSlider.value
        numIterations = iterationsSlider.value
        stoppingCriterion = stoppingCriterionSlider.value
        ourPerceptron.setWs([startingWeight.value], startingBias.value)
        output = gradientDescent(ourPerceptron, mseCost, np.reshape(x[0:50], (-1, 1)),
                                 y[0:50], learningRate, numIterations, stoppingCriterion)
        display(compareOutputs(x, y, ourPerceptron, xRange))
        display(plotOutputs(output))

# Set the function to be called when the button is clicked
updateBtn.on_click(updateClick)

# Display the widgets
display(learningRateSlider, iterationsSlider, stoppingCriterionSlider,
        startingWeight, startingBias, updateBtn, outputWidget)
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...

Résumé

  • Un perceptron est une somme pondérée plus un biais, passée dans une fonction d’activation.
  • La règle d’apprentissage classique du perceptron ne converge que sur des données linéairement séparables. Elle a séparé nos fenêtres d’événements et de bruit à partir de deux caractéristiques de détecteur, et ce choix de deux caractéristiques nous a permis de tracer la frontière de décision.
  • Sans activation, le même perceptron est un régresseur linéaire. La descente de gradient sur la perte MSE retrouve presque la même droite que les OLS.
  • Le taux d’apprentissage règle le compromis entre convergence lente et divergence.

La suite : 4.1 Réseaux de neurones empile des perceptrons en couches pour traiter des données qu’une seule droite ne peut pas séparer.