Python Notebook Series 1

28 augustus, 2020

Rick Vink

“Een serie geschreven in Notebook stijl. Speciaal voor lezers die algoritmes graag willen toepassen op eigen datasets.”

Deel 1: Vier manieren om AI algortimes te optimaliseren.

Wanneer je een Machine Learning/AI algoritme wilt maken zijn er meerdere beslissingen die je moet maken. Welk doel geef ik mijn AI? Hoe voed ik het algoritme data? Welke data? Welk model is geschikt? Hoe gaat het model zichzelf optimaliseren? Voor elk van deze aspecten zijn er vele variaties bedacht. Voor de laatste vraag zijn er hele interessante technieken bedacht die gebasseerd zijn op natuurlijke evolutie en de manier waarop het brein leert. Tijdens deze editie van de Python Notebook series kijken we naar vier verschillende strategieën voor het optimaliseren van het model:

Grid Search
Normal Equation
Gradient Decent
Genetic Algorithm (evolutie)

Om de kracht van iedere strategie te demonstreren, zullen we de verschillende strategieën gebruiken voor het oplossen van een specifiek probleem.

Probleemstelling

Stel dat er een trein is waarvan we de versnelling willen bepalen wanneer we een beperkt aantal metingen hebben als data. Specifieker: stel we hebben 1000 metingen gedaan over een tijdsduur van 5 seconden, waarbij de metingen onderhevig zijn aan meetonzekerheden ten gevolge van ruisbronnen.

In [1]:

#Allereerst beginnen we met het importeren van de benodigde libraries

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import random
%matplotlib notebook

Genereer dataset

Om de dataset te genereren maken we gebruik van de NumPy’s random library. We genereren data vanuit een meetinstrument dat 1000 metingen heeft gedaan in een tijdspanne van 5 seconden. De eigenschappen van ons meting zijn als volgt:

Er is een error doordat de meting niet precies bij het begin van de acceleratie van de trein startte (offset error).
Er is een ruis door verschillende ruisbronnen die invloed hadden op de meting.
Het meetinstrument meet niet in gelijke tijdstappen.

In de code hierop volgend simuleren we de data vanuit de hierboven beschreven meting.

In [2]:

number_of_points = 1000     # aantal meetpunten
total_time = 5              # Totale tijd waarover we metingen doen
offset = 2                    # Een offset error die we optellen (te laat begonnen met meten)
acceleration = 3            # De echte waarde van de versnelling van de trein
noise = 0.5                 # Er is een ruisbron aanwezig met waarde 0.5 
# Maak 1000 willekeurige tijd punten aan: het meetinstrument meet niet in gelijke tijdstappen. 
timepoints = np.random.random(number_of_points)*total_time
# Soorteer de tijd punten van laag naar hoog
timepoints = np.sort(timepoints)
# Creer ruis per meting van de grootte willekeurig per meting
noise = noise * np.random.randn(number_of_points)
# Maak een lijst van de zogenaamde gemeten waardes
speeds = acceleration * timepoints + noise + offset

Laten we de data visualiseren. Er is goed te zien dat de data niet mooi op een rechte lijn ligt, vanwege de ruisbron.

In [3]:

figure = plt.figure()
plt.scatter(timepoints, speeds)
plt.ylabel('Snelheid (m/s)')
plt.xlabel('Tijd (s)')
plt.show()

Model kiezen

We zien al dat het een rechte lijn is dus laten we uit gaan van het volgende model:

v=a⋅t+v0

Waarbij v de snelheid is, a de versnelling, t de tijd, en v_0 de beginsnelheid toen de meting begon. De snelheid en tijd weten we uit onze metingen. De versnelling en de beginsnelheid zijn waarden die we met ons model willen achterhalen.

Het doel

Het doel is om een zo goed mogelijke schatting te maken voor de versnelling. Maar hoe moeten we dit meten? Stel we willen weten wat de fout is van de waarde a = 1 en v_0 = 1.

In [4]:

# Onze voorspelde waarde voor de versnelling (die dus duidelijk fout is)
acceleration_prediction = 1
bias_prediction = 1
# Voorspelde snelheden
speed_predictions = acceleration_prediction * timepoints + bias_prediction

In [5]:

plt.plot(timepoints, speed_predictions, color='red')
plt.scatter(timepoints, speeds)
plt.ylabel('Snelheid (m/s)')
plt.xlabel('Tijd (s)')
plt.legend(['Voorspelling','Werkelijkheid'])

Out[5]:

We zien dat we erg ver af zitten van de werkelijkheid. Graag willen we deze fout kwantificeren in een getal zodat we dat kunnen vergelijken met andere schattingen. Dit doen we door elke keer het verschil van onze voorspelling en de werkelijke waarde te nemen en daar het kwadraat van te nemen. Dit is een veel gebruikte methode om te meten hoe ver we af zitten van de werkelijke waarde.

F_out_i = (v_verschil_i)^2 = (v_voorspelling_i - v_echt_i)^2 F_out = \sumF_out_i

In [6]:

# Laten we het verschil bepalen van een van de meetpunten 
sample_number = 500

# laten we het tijdpunt, de echte snelheid en de voorspelde snelheid pakken
time_point = timepoints[sample_number]
speed_true = speeds[sample_number]
speed_predicted = speed_predictions[sample_number]

# maake een figuur
plt.plot(timepoints, speed_predictions, color='red')
plt.plot([time_point, time_point], [speed_true, speed_predicted],c='orange')
plt.scatter(timepoints, speeds)

plt.ylabel('Snelheid (m/s)')
plt.xlabel('Tijd (s)')
plt.legend(['Voorspelling','Verschil','Werkelijkheid'])

Out[6]:

Vinden van de optimale waarde

Nu zijn we bij het punt dat we moeten bepalen hoe we onze waarde van de versnelling willen gaan aanpassen zodat we uiteindelijk in de buurt komen van de echte waarde. We gaan hiervoor vier methodes proberen: Grid Search, Normal Equation, Gradient Descent en Genetic Algorithm.

Grid Search

Een oplossing kan zijn om vele verschillende waardes van de versnelling te proberen en kijken welke de kleinste fout heeft. Dit klinkt heel eenvoudig en lomp (en dat is het ook) maar het maakt wel een mooi figuurtje.

In [7]:

# De range aan vernellingen die we gaan testen
lowest_testing_acceleration = 0
highest_testing_acceleration = 10
# De range aan beginsnelheden die we gaan testen
lowest_bias = -5
highest_bias = 5
# We testen meerdere waardes tussen de uiterste
number_of_accelerations_to_test = 300
number_of_biases_to_test = 200
# De verschillende versnellingen die we gaan proberen
acceleration_tests = np.linspace(lowest_testing_acceleration,highest_testing_acceleration,number_of_accelerations_to_test)
bias_tests = np.linspace(lowest_bias,highest_bias,number_of_biases_to_test)

In [8]:

def calculate_loss_speed(acceleration_tests, bias_tests):
    # We gaan voor elk tijdpunt berekenen wat de fout is tussen de voorspelling
    # en de werkelijke waarde. Dit doen we voor elke combinatie van de versnelling
    # en de begin snelheid. Daarom maken we een nieuwe lege matrix ter grote van
    # elke combinatie van de acceleratie en de begin snelheid
    sum_losses = np.zeros(acceleration_tests.shape)
    for timepoint_index, timepoint in enumerate(timepoints):
        # Bereken de fout op dat tijdpunt voor elke combinatie van de versnelling
        # en de begin snelheid
        speed_predict = acceleration_tests * timepoint + bias_tests
        # Tel het op bij de huidige stand
        sum_losses += (speed_predict - speeds[timepoint_index])**2
    # Bereken de gemiddelde loss omdat dat wat netter is
    mean_losses=sum_losses/len(timepoints)
    return mean_losses

# Maak een grid van elke combinatie van de versnelling en de begin snelheid
A, B = np.meshgrid(acceleration_tests, bias_tests)
# De berekende loss voor elke combinatie van de versnelling en de begin snelheid 
L = calculate_loss_speed(A, B)

In [10]:

# Laten we kijken naar het resultaat
figure = plt.figure(figsize=(14,14))
ax = plt.axes(projection='3d')
ax.plot_surface(A, B, L, rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')
ax.set_title('surface');
ax.set_ylabel('begin snelheid (m/s)');
ax.set_xlabel('Versnelling (m/s^2)');
ax.set_zlabel('Loss');

ax.plot([lowest_testing_acceleration-1 ,highest_testing_acceleration+1],[offset,offset],[0,0], linewidth=10)
ax.plot([acceleration,acceleration],[lowest_bias-1 ,highest_bias+1],[0,0], linewidth=10)

Out[10]:

In [11]:

# Vind de minimale waarde van de versnelling
min_accelerations = np.min(L, axis=0)
min_acc_arg = (np.argmin(min_accelerations))
min_acc = acceleration_tests[min_acc_arg]
print('Versnelling van laagste loss: {}'.format(min_acc))

# Vind de minimale waarde van de bias
min_biases = np.min(L, axis=1)
min_bias_arg = (np.argmin(min_biases))
min_bias = bias_tests[min_bias_arg]
print('Beginsnelheid van laagste loss: {}'.format(min_bias))

Versnelling van laagste loss: 3.0100334448160533
Beginsnelheid van laagste loss: 1.9849246231155782

Zoals je ziet komen de voorspelde waardes zeer in de buurt van wat we hadden verwacht.

In [12]:

acceleration_prediction = min_acc
bias_prediction = min_bias
# Voorspelde snelheden
speed_predictions = acceleration_prediction * timepoints + bias_prediction

plt.plot(timepoints, speed_predictions, color='red')

plt.scatter(timepoints, speeds)
plt.ylabel('Snelheid (m/s)')
plt.xlabel('Tijd (s)')
plt.legend(['Voorspelling','Werkelijkheid'])

Out[12]:

Met deze methode halen we een zeer goede schatting van de versnelling van de trein.

Normal equation

Een andere zeer logische oplossing voor dit probleem is om onze wiskunde/natuurkunde kennis direct toe te passen.

v=a\cdott+v_0 a=(v-v_0)/t

Deze methode werkt echter alleen in gevallen waarin het probleem analytisch op te lossen is. Dit maakt gebruikt met de zogeheette normal equation waarmee je direct het antwoord kan uitrekenen:

variables=(X^T X)^-1 X^T y

Waarbij X de input parameter matrix is en de y de vector met de resultaten. In ons geval hebben we eigenlijk maar een input parameter en ook nog eens een bias die we hierin ook moeten verwerken. Dit doen we door een constane waard evan 1 ook toe te voegen bij onze input parameters.

In [13]:

# Onze y
y =  speeds
# Onze X Waarbij we een constante waarde van 1 er bij doen om zo de bias mee te nemen
X = np.concatenate((
        np.ones((timepoints.shape[0],1)), 
        timepoints.reshape(timepoints.shape[0],1))
    ,axis=1)

In [14]:

# Zo zien de eerste 5 waardes er uit van de y term
y[0:5]

Out[14]:

array([3.04644207, 1.91328379, 1.38991765, 1.98514308, 1.5779907 ])

In [15]:

# En zo zien de eerste 5 waardes er uit van onze X term
X[0:5,:]

Out[15]:

array([[1.        , 0.0037286 ],
       [1.        , 0.00696756],
       [1.        , 0.01468777],
       [1.        , 0.01776416],
       [1.        , 0.02093316]])

In [16]:

# Nu passen we de normal equation toe
variables = np.dot(np.linalg.inv(np.dot(X.T, X)), np.dot(X.T, y))
variables

Out[16]:

array([2.00286255, 3.00240928])

Het antwoord hierboven ligt zeer dicht bij gezochte waarden. De methode van de Normal Equation is toe te passen op problemen die met een simpele formule te beschrijven zijn, maar niet voor complexere problemen.

Gradient Descent

We zagen bij de Grid Search methode dat er een hele mooie helling is die leidt naar de minimale fout. Het is dan ook een zeer populaire methode om de helling te bepalen (de zogehete gradient) en elke keer een stapje te doen naar beneden. We nemen steeds een kleine stap de helling af dat geschaald is met een learning rate. Hoe hoger de learning rate, hoe sneller we naar beneden gaan. Kiezen we echter de learning rate te hoog, dan is mogelijk dat we dan juist voorbij de minimale waarde schieten.

anew=aprevious-learning_rate\cdotdFoutda

v0,new=v0,previous-learning_rate\cdotdFoutdv0

Voor deze methode nemen we een willekeurige waarde voor onze versnelling en berekenen we de helling die daarbij hoort.

dFoutda=d(vvoorspelling-vecht)2da=d(a\cdott+v0-vecht)2da=2(a\cdott+v0-vecht)\cdott=2(vvoorspelling-vecht)\cdott

dFoutdv0=d(vvoorspelling-vecht)2dv0=d(a\cdott+v0-vecht)2dv0=2(a\cdott+v0+vecht)=2(vvoorspelling-vecht)

Dit process wordt iteratief uitgevoerd waarbij we beginnen met een voorspelling van de waarde van de versnelling en steeds een aanpassing maken totdat we in de buurt komen van een goede voorspelling.

In [17]:

# De leersnelheid bepaald hoe groot onze stapgrootte is
learning_rate = 0.01
# Begin punt van de versnelling en de begin snelheid
acceleration_prediction = 10
bias_prediction = 5
# Het aantal itteraties die we willen doen
number_of_itteration = 100

# Laten we de resultaten bij houden
gd_mean_losses = []
gd_a = []
gd_b = []

for step in range(number_of_itteration):
    # Bereken voorspelde snelheden v_voorspeld = a * t
    speed_predictions = acceleration_prediction * timepoints + bias_prediction
    # Bereken de gradient: (2 v_voorspeld - v_echt) * t
    slope_a = np.mean(2 * (speed_predictions - speeds) * timepoints)
    # Bereken de gradient: (2 v_voorspeld - v_echt)
    slope_b = np.mean(2 * (speed_predictions - speeds))
    # Laten we ook maar de fout berekening (optioneel)
    mean_loss = np.mean((speed_predictions - speeds)**2)
    gd_mean_losses.append(mean_loss)
    gd_a.append(acceleration_prediction)
    gd_b.append(bias_prediction)
    # Pass de voorspelling van de versnelling aan
    acceleration_prediction = acceleration_prediction - learning_rate * slope_a
    bias_prediction = bias_prediction - learning_rate * slope_b

In [18]:

plt.plot(gd_mean_losses)

Out[18]:

In [19]:

# Laten we kijken naar het resultaat
figure = plt.figure(figsize=(14,14))
ax = plt.axes(projection='3d')
ax.plot_surface(A, B, L, rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')
ax.set_title('surface');
ax.set_ylabel('begin snelheid (m/s)');
ax.set_xlabel('Versnelling (m/s^2)');
ax.set_zlabel('Loss');

ax.plot(gd_a,gd_b,gd_mean_losses, linewidth=10)

Out[19]:

En de fout bij iedere iteratie

In [20]:

plt.scatter(range(len(gd_mean_losses)), gd_mean_losses)
plt.ylabel('Fout')
plt.xlabel('Itteraties')

Out[20]:

En zo ook de voorspelde waarde van de versnelling

In [21]:

plt.scatter(range(len(gd_a)), gd_a)
plt.scatter(range(len(gd_b)), gd_b)
plt.legend(['a (m/s^2)','beginsnelheid (m/s)'])
plt.xlabel('Itteraties')

Out[21]:

Zoals je ziet heb je hiermee na een paar stappen al een zeer goede voorspelling van de versnelling en beginsnelheid.

Genetic Algorithm

Last but not least, Genetic Algorithms. De naam zegt het al, dit algorithme is gebaseerd op natuurlijke selectie wat een drijvende kracht is van biologische evolutie. Het grote verschil tussen genetic algorithme en Gradient Descent is dat we in dit geval meerdere versies van de code hebben in parallel en na een stap bepalen welke versies we selecteren en welke afvallen.

In [22]:

# We gaan 100 versies van de code naast elkaar laten leven
number_of_codelets = 100 # codelet: schattige naam voor een versie van de code in de zelfde generatie.
# De versies van de code gaan we voor 20 generaties volgen
number_of_generations = 20

# Laten we verschillende codelets maken waarbij elke een andere waarde voor de versnelling en beginsnelheid is.
# Elke codelet is een waarde voor de versnelling van voornamelijk tussen de -10 en 10 voor de snelheden en -5 en 5 voor de beginsnelheden.
# De waardes kunnen hier ook buiten liggen aangezien het een Gauss distributie is
begin_accs  = (20 * np.random.rand(number_of_codelets) - 10) # versnellings termen
begin_biases = (10 * np.random.rand(number_of_codelets) - 5)  # beginsnelheid termen
codelets = []
for begin_acc, begin_bias in zip(begin_accs, begin_biases):
    codelets.append(np.array([begin_acc, begin_bias]))

# We gaan steeds de waarde bij houden hoe goed de codelets het doen
fitness_values_results = []
# En we onthouden steeds de waarde van elke codelet
results = []

# Meerdere generaties codelets
for generation in range(number_of_generations):
    # Elke codelet zal moeten worden getest
    # Waarbij elke codelet met een waare van 0 begint (die we laten gaan overschrijven)
    fitness_values = np.zeros(number_of_codelets)

    # We gaan elke codelet testen
    for codelet_index, codelet in enumerate(codelets):
        
        # Bereken alle snelheden voor elk tijdspunt volgens een codelet
        speed_predictions = codelet[0] * timepoints + codelet[1]
        
        # Fitness testen
        # Hoe lager de loss hoe fitter de functie
        loss = np.mean((speed_predictions - speeds)**2)
        fitness_value = -loss
        
        fitness_values[codelet_index] = fitness_value
#         print(fitness_value)
        
    # Sla de waardes op voor analyse
    fitness_values_results.append(fitness_values)
    results.append(codelets)
        
    # Natuurlijke selectie
    # We gaan semi willekeurig een paar codelets laten leven.
    # Hierbij geeft een hogere fitness waarde een hogere kans van evolueren
    survival_probability = fitness_values - np.min(fitness_values)
    survival_probability /= np.sum(survival_probability)
    survivors_indices = np.random.choice(number_of_codelets, int(number_of_codelets/2), p=survival_probability, replace=False)
    # Voeg de overlevende toe aan de nieuwe generatie
    survivors = []
    for survicor_index in survivors_indices:
        survivors.append(codelets[survicor_index])
    
    # Maak gemuteerde clonen van de overlevende
    mutated_babies = survivors + 0.1 * np.random.randn(number_of_codelets - int(number_of_codelets/2),2)
    
    codelets = np.concatenate((survivors, mutated_babies))

Laten we kijken naar het resultaat van het algorithme. Hierbij plotten we bij iedere generatie steeds de beste, gemiddelde en slechtste waarde.

In [24]:

# Bepaal de beste van elke generatie
best_nodelet = np.max(fitness_values_results,axis=1)
# Bepaal degene in het midden van elke generatie
mean_nodelet = np.median(fitness_values_results,axis=1)
# Bepaal de slechtste van elke generatie
worst_nodelet = np.min(fitness_values_results,axis=1)
plt.plot(best_nodelet)
plt.plot(mean_nodelet)
plt.plot(worst_nodelet)
plt.ylabel('Fout')
plt.xlabel('Generatie')
plt.legend(['Beste','Gemiddelde','Slechtste'])

Out[24]:

En zo zien we ook dat de best performende waarde ook vrij snel richting het juiste antwoord komt.

In [25]:

results = np.array(results)
best_nodelet_indeces = np.argmax(fitness_values_results,axis=1)
best_accelerations = results[range(number_of_generations),best_nodelet_indeces]
plt.plot(best_accelerations)
plt.legend(['Versnelling','beginsnelheid'])
plt.xlabel('Generatie')

Out[25]:

Welke moet ik kiezen?

De vier methodes zijn allemaal vrij divers en hebben allemaal zo hun voor en nadelen.

Grid Search

De gridsearch valt al heel snel af wanneer het probleem al wat meer complex wordt. Het is namelijk vrij computational expensive en het vooraf bepalen van de breedte en hoeveelheid samples van het grid is niet altijd triviaal.

Analytisch

Deze methode is een zeer logische optie wanneer het model analytisch oplosbaar is.

Gradient Descent

Deze methode is zeer populair. Er zijn veel voordelen voor het gebruiken van deze methode:

Het is vrij effecient omdat het de gradient nauw volgt.
Is toepasbaar op complexe systemen.

Echter zijn er ook nadelen aan Gradient Descent:

De meest voordehandliggende vereiste is dat er een gradient berekent moet kunnen worden bij ieder punt. Dit is echter niet altijd het geval (zoals bijvoorbeeld bij een stap functie).
Aanpassingen aan de structuur van het model kunnen niet worden meegenomen tijdens het leren. Hier kan je denken aan het breder of dieper maken van een neuraal netwerk.

Genetic Algortithm

Genetische Algorithms zijn wat minder populair geworden de laatste jaren met de komst van neurale netwerken. Dit omdat Gradient Descent hier zeer goed op werkt en efficienter is. Dit betekent echter niet dat het geen plek meer heeft in onze toolbox. Er zijn nog genoeg voordelen om genetische algorithms toe te passen:

Hyper parameter tuning. Het maakt een genetisch algoritme niet uit wat je aanpast. Dit kunnen trainbare variabelen zijn zoals we in dit artikel deden, maar het kunnen ook hele aanpassingen zijn aan je model of het aanpassen van learning rates etc.
Goed te parallelliseren. Genetische algoritmes kunnen prima naast elkaar werken op andere computers gevolgd met een selectie proces.

Combinatie van Gradient Descent en Genetic Algorithm

De ene methode sluit het andere niet uit. Het combineren van Gradient Descent en Genetic Algorithm is een hele sterke samenwerking. Zo kan je elke generatie meerdere computers een neuraal netwerk laten trainen gevolgd met een fitness test van elk model, waarna je de modelstructuren kan aanpassen en weer een nieuwe generatie kan laten trainen. Het grote voordeel hiervan is dat je dan de hyper parameters ook kan mee nemen in je trainings process van je netwerk.

Conclusie

Zoals het vaak zo is heeft elke methode wel weer zo zijn toepassing. Voor visualisatie is Grid Search leuk, kleine modelen kan je prima doorrekenen met de Normal Equation en voor complexe neurale netwerken kan je Gradient Descent toepassen met een vleugje van Genetic Algorithm hier en daar.

In dit artikel heb ik één vorm van Gradient Descent en van Genetic Algorithm laten zien. Hier zijn ook allemaal varianten van waar over ook op zichzelf een artikel over geschreven kan worden.