Es stehen mittlerweile viele Werkzeuge wie Keras, als Applikationsschnittstelle für Bibliotheken wie tensorflow, zur Verfügung um neuronale Netzwerke zu programmieren. Um ein besseres Verständnis der Zusammenhänge zu erlangen, möchte ich jedoch ein neuronales Netzwerk von Grund auf, d.h. im Wesentlichen nur unter Verwendung von Phyton und numpy programmieren.
Ein neuronales Netzwerk, genauer gesagt ein Feedforward Neural Network (FNN), besteht aus einer Eingabeschicht, einer oder mehreren verborgener Schichten und einer Ausgabeschicht.
Die einzelnen Schichten \(A\) beinhalten sogenannte Neuronen, die über lineare Abbildungen
$$ Z^{l+1}=W^{l+1}\cdot A^l + B^{l+1}$$
miteinander verknüpft sind. Dabei bezeichnet \(W^{l+1}\) eine Matrix, deren Dimension der Zahl der Neuronen der Schicht \(l+1\) mal der Zahl der Neuronen der Schicht \(l\) entspricht. \(B^{l+1}\) bezeichnet einen Vektor der addiert wird. Sollen Bilder verarbeitet werden, stellen die Neuronen z.B. die einzelnen Pixel des Bildes dar. Die Matrizen \(W\) und \(B\) werden im Rahmen eines Trainings so angepasst, dass die Ausgabeschicht das gewünschte Ergebnis möglichst gut erreicht. Da durch lineare Abbildungen nur lineare Zusammenhänge 'angefittet' werden können, werden die errechneten \(Z \) Werte jeder Schicht durch nicht lineare Funktionen 'aktiviert', und dadurch in Werte von \(A\) abgebildet. Die einfachste, und sehr häufig verwende Aktivierungsfunktionen ist die sogenannte Rectified Linear Unit (ReLU) Funktion. Die ReLu Funktion liefert 0 wenn \(z<0\) und \(z\) wenn \(z\geq 0\).
Sollen z.B. Bilder klassifiziert werden, dient diese Funktion im Netzwerk dazu die Ausgabe in eine Wahrscheinlichkeitsverteilung über die vorherzusagenden Ausgabeklassen umzuwandeln. Für Mehrklassenprobleme wird im Allgemeinen die Softmax-Funktion verwendet.
\[ \text{softmax}(z_j)= s_j = \frac{e^{z_j}}{\sum_{i=1}^{k} e^{z_i}} \]
Dabei ist \(z\) der Eingabevektor mit \(k\) Komponenten, \(z_j\) die \(j\)-te Komponente von \(z\), und die Summe im Nenner läuft über alle \(k\) Komponenten. Grafisch lässt sich die Funktion in Abhängigkeit einer Komponente wie folgt darstellen:
Die Funktion kann als normalisierte Exponentialfunktion betrachtet werden, und ist eine Verallgemeinerung der logistischen Funktion auf mehrere Dimensionen. Bei welchem \(z\) Wert der 'Übergang' von 0 nach 1 einer Komponente stattfindet, hängt von den Werten der anderen Komponenten ab. Durch die Summe im Nenner ist jedoch stets gewährleistet das die Summe aller Einzelwahrscheinlichkeiten 1 beträgt.
Das Training des neuronalen Netzes läuft nun wie folgt ab. Zu Beginn des ersten Durchlaufs werden die Matrizen \(W\) und \(B\) zufällig initialisiert. Anschließend erfolgt die sogenannte 'forward propagation' mit aufbereiteten Trainingsdaten. Die errechneten Wahrscheinlichkeiten erreichen in der Regel nicht 100% für eine Kategorie. Wie stark die Ergebnisse abweichen, kommt in der Kostenfunktion zum Ausdruck. Soll mehr als eine Kategorie bestimmt werden, empfiehlt es sich die sogenannte Kreuzentropie zweier Wahrscheinlichkeitsvektoren als Kostenfunktion zu verwenden.
$$ L(y,s) = - \sum_{i=1}^{C} y_i \cdot \log(s_i) $$
Dabei bezeichnet \( C \) die Anzahl der Klassen, \(s_i\) die Ausgabewerte der letzten Schicht, und \(y_i\) den Ergebnisvektor in 'one hot' oder '1 aus n' Kodierung. Je kleiner die Kreuzentropie, desto ähnlicher sind die beiden Wahrscheinlichkeitsverteilungen, denn \(-log(s_i)\) wird umso kleiner je mehr sich \(s_i\) der 1 annähert.
Im weiteren Verlauf des Trainings werden die Parameter \(W\) und \(B\) nach der sogenannten 'gradient descent' Methode optimiert. Dazu wird die Änderung der Kostenfunktion für jeden Layer und jeden Parameter bestimmt. Anschließend wird das neuronale Netzwerk in der 'back propagation' durchlaufen, und die Parameter für den nächsten Trainingsdurchlauf angepasst. Der letzte nicht aktivierte Layer \(ll\) soll sich also iterativ wie folgt ändern:
$$ Z^{ll}:=\nabla_{Z^{ll}}L = Z^{ll}-\alpha\frac{dL}{dZ^{ll}}$$
Dabei bezeichnet \(\alpha\) die sogenannte Lernrate. Um die Ableitung der Kostenfunktion \(L\) nach \(Z\) bestimmen zu können, müssen wir zunächst die Softmaxfunktion ableiten. Das erfordert etwas Aufwand, aber wenn wir die partielle Ableitung des Logarithmus bilden, können wir die Rechnung vereinfachen.
$$ \frac{\partial}{\partial z_j} \log(s_i) = \frac{1}{s_i} \cdot \frac{\partial s_i}{\partial z_j}, $$
Der Ausdruck auf der rechten Seite folgt direkt aus der Kettenregel. Als Nächstes ordnen wir die obere Formel um.
$$\frac{\partial s_i}{\partial z_j} = s_i \cdot \frac{\partial}{\partial z_j} \log(s_i)$$
Die linke Seite ist genau die partielle Ableitung, die wir suchen. Wir betrachten zuerst den Logarithmus von s.
$$\begin{align*}\log s_i & = \log \left( e^{z_i} / \sum_{i=1}^n e^{z_i} \right) = z_i - \log \left( \sum_{i=1}^n e^{z_i} \right)\end{align*}$$
Die partielle Ableitung des resultierenden Ausdrucks ist:
$$\frac{\partial}{\partial z_j}\log s_i = \frac{\partial z_i}{\partial z_j}-\frac{\partial}{\partial z_j}\log \left( \sum^n_{i=1}e^{z_i}\right)$$
Der erste Term auf der rechten Seite
$$ \frac{\partial z_j}{\partial z_j}= \begin{cases}1, & i = j \\ 0, & \text{sonst} \end{cases}$$
kann prägnant mit der Indikatorfunktion \(1\{i=j\}\) geschrieben werden. Die Indikatorfunktion nimmt den Wert 1 an, wenn \(i=j\) ist, und sonst 0. Der zweite Term auf der rechten Seite kann durch Anwendung der Kettenregel ausgewertet werden.
$$\frac{\partial}{\partial z_j} \log s_i = 1\{i = j\} - \frac{1}{\sum_{l=1}^n e^{z_l}} \left( \frac{\partial}{\partial z_j} \sum_{l=1}^n e^{z_l}\right)$$
Die partielle Ableitung der Summe zu erhalten, ist trivial.
$$\frac{\partial}{\partial z_j} \sum_{i=1}^n e^{z_i} = \frac{\partial}{\partial z_j} [e^{z_1} + e^{z_2} + ... + e^{z_j} + ... + e^{z_n}] = \frac{\partial}{\partial z_j} [e^{z_j}] = e^{z_j}$$
Eingesetzt in die obige Formel führt dies zu:
$$\frac{\partial}{\partial z_j} \log s_s = 1\{i=j\} - \frac{e^{z_i}}{\sum_{i=1}^n e^{z_i}} = 1\{i=j\}-s_j$$
Schließlich müssen wir den obigen Ausdruck mit \(s_i\) multiplizieren wie am Anfang des Abschnitts gezeigt:
$$\frac{\delta s_i}{\delta z_j}=s_i \frac{\partial}{\partial z_j}\log(s_i) =s_i \left( 1\{i=j\}-s_j \right)$$
Damit ist unsere Ableitung abgeschlossen. Wir haben eine Formel für alle Elemente der Jacobi-Matrix \(J\) erhalten. Für den Spezialfall von n = 4 gilt zB.:
$$J_{\text{softmax}} = \begin{bmatrix}s_1 \cdot (1 - s_1) & -s_1 \cdot s_2 & -s_1 \cdot s_3 & -s_1 \cdot s_4 \\ -s_2 \cdot s_1 & s_2 \cdot (1 - s_2) & -s_2 \cdot s3 & -s_2 \cdot s_4 \\ -s_3 \cdot s_1 & -s_3 \cdot s_2 & s_3 \cdot (1 - s_3) & -s_3 \cdot s4 \\ -s_4 \cdot s_1 & -s_4 \cdot s_2 & -s_4 \cdot s_3 & s_4 \cdot (1 - s_4) \end{bmatrix}$$
Als nächstes wollen wir nun die Ableitung der Kostenfunktion in Bezug auf die Eingabe \(z\) der letzten Aktivierungsschicht berechnen.
$$\frac{\partial L}{\partial z_j} = - \frac{\partial}{\partial z_j}\sum_{i=1}^{c} y_i \cdot \log(s_i) = - \sum_{i=1}^{c} y_i \cdot \frac{\partial}{\partial z_j}\log(s_i) = -\sum_{i=1}^{c} \frac{y_i}{s_i} \cdot \frac{\partial s_i}{\partial z_j}$$
Fügen wir die Ableitungen ein, die wir oben erhalten haben,
$$\frac{\partial L}{\partial z_j} = -\sum_{i=1}^{c} \frac{y_i}{s_i} \cdot s_i \cdot (1\{i=j\}-s_j) = -\sum_{i=1}^{c} y_i \cdot (1\{i=j\}-s_j)$$
und multiplizieren aus
$$\frac{\partial L}{\partial z_j}= \sum_{i=1 }^c y_i \cdot s_j-\sum_{i=1 }^c y_i \cdot 1\{i=j\})$$
dann ergibt sich
$$\frac{\partial L}{\partial z_j} = \sum_{i=1}^{c} y_i \cdot s_j - y_j$$
Als nächstes ziehen wir \(s_i\) aus der Summe, da es nicht vom Index \(i\) abhängt
$$\frac{\partial L}{\partial z_j} = s_j \cdot \sum_{i=1}^c y_i -y_j = s_j - y_j$$.
In prägnanter Vektorschreibweise erhalten wir:
$$\nabla_{Z^{ll}}L = S -Y$$
Als nächstes wollen wir bestimmen, wie sich \(W\) nach einem Trainingsdurchlauf ändern soll. Für die Änderung von \(W\) soll analog zur Änderung von \(Z\) gelten:
$$W^l:=W^l-\alpha '\frac{dL}{dW^l}=W^l-\alpha '\frac{dL}{dZ^l}\cdot\frac{dZ^l}{dW^l}$$
Der Ausdruck \(dZ/dW\) lässt sich mathematisch ähnlich schwer formulieren wie der bereits berechnete Ausdruck \(dL/dZ\). Es gilt:
$$\frac{\partial z_{ik}^l}{\partial w_{mn}^l}=\frac{\partial}{\partial w_{mn}^l}\sum_{j=1}^{\text{dim(layer l)}}w_{ij}^l a_{jk}^{l-1} + b_{ik}^l=1\{m=i,n=j\}a_{nk}^{l-1}$$
Dabei bezeichnet \(a_{nk}\) das Element einer Matrix der Dimension \(dim(\text{layer l}) * N\), wobei \(N\) die Zahl der Trainingsdatensätze darstellt. Durch die Verwendung des Skalarprodukt \(dL/dZ \cdot dZ/dW\) möchten wir die Änderung in \(W\) über alle Datensätze mitteln. Wir normieren daher auf \(1/N\) und transponieren \(A^{l-1}\). Es gilt:
$$\lambda ' \left( \frac{dL}{dZ^l}\cdot \frac{dZ^l}{dW^l} \right)_{mn}=\frac{\lambda}{N}\sum_{k=1}^N \left( \frac{dL}{dz^l}\right)_{mk}\cdot a_{kn}^{l-1}$$
oder übersichtlicher
$$\nabla_{W^l}L = \frac{1}{N} \nabla_{Z^l}L\cdot {A^{l-1}}^T$$
Auch für die iterative Änderung von \(B\) gilt:
$$B^l:=B^l-\alpha '\nabla_{B^l}L=B^l-\alpha '\frac{dL}{dB^l}= B^l-\alpha '\frac{dL}{dZ^l}\cdot\frac{dZ^l}{dB^l}=B^l-\alpha '\frac{d L}{dZ^l}^l\cdot 1$$
Dabei bezeichnet \(1\) eine Einheitsmatrix der Dimension \(dim(\text{layer l}) * N\). In Vektorschreibweise gilt wieder:
$$\nabla_{B^l}L=\frac{1}{N}\nabla_{Z^l}L\cdot 1^T$$
Wir multiplizieren erneut mit der transponierten (Einheits)matrix um über alle Datensätze zu mitteln.
Zu guter Letzt gilt für das update der inneren Schichten \(Z\):
$$Z^l:=Z^l-\alpha '\nabla_{Z^l}L = Z^l-\alpha '\frac{dL}{dZ^l}=Z^l-\alpha '\frac{dZ^{l+1}}{dA^l}\cdot\frac{dL}{dZ^{l+1}}\cdot\frac{dA^l}{dZ^l}$$
und in Vektorschreibweise bedeutet das:
$$\nabla_{Z^l}L={W^{l+1}}^T\cdot\nabla_{Z^{l+1}}L\odot\text{ReLu}'(Z^l)$$
Dabei steht \(\odot\) für die elementweise Multiplikation bzw. das Hadamard Produkt und \(ReLu '\) für die Ableitung der ReLu Aktivierungsfunktion. Diese Ableitung ist trivialerweise durch
$$ReLu'(z) = \begin{cases} 1 & \text{wenn } z > 0 \\0 & \text{wenn } z \leq 0 \end{cases}$$
gegeben. Damit sind die mathematischen Grundlagen unseres neuronalen Netzwerkes im Wesentlichen erarbeitet, und wir können uns an die Programmierung wagen.
Unser neuronales Netzwerk soll als Trainingsdaten den berühmten MNIST Datensatz verwenden. Dieser besteht aus 42000 handgeschriebenen Ziffern in 8 bit Graustufen mit der Auflösung 28 x 28 Pixel. Die Eingabeschicht \(A^0\) besteht daher aus 784 Einheiten, die den 784 Pixeln in jedem 28 x 28-Eingabebild entsprechen. Der Programmcode ist so gestaltet, dass beliebig viele innere Schichten eingefügt werden können. Bereits mit nur einer inneren Schicht lassen sich jedoch sehr gute Ergebnisse erzielen. Diese eine versteckte Schicht \(A^1\) besteht aus 10 Einheiten mit ReLU-Aktivierung, und die Ausgabeschicht \(A^2\) besteht ebenfalls aus 10 Einheiten, die den zehn Ziffernklassen mit Softmax-Aktivierung entsprechen.
Der gesamte Code steht als Kaggle notebook zum Ausprobieren zur Verfügung. Ich habe versucht die einzelnen Schritte möglichst gut im Programmcode zu kommentieren.
import numpy as np import pandas as pd import time from datetime import datetime from matplotlib import pyplot as plt # set structure of neural network layer_dim = [784,10,10] # number of neurons per layer (including input and output layer) L = len(layer_dim) # number of layers LL = L-1 # index of Last Layer # set training parameter test_size = 1000 # set number of test samples, note: test dataset will be removed from training dataset training_rate = 0.10 # set training rate iterations = 200 # set number of iterations # set global variables params = {} # create dictionary for parameter W and b cache = {} # create dictionary for A and Z grads = {} # create dictionary for parameter gradients cost = np.zeros(iterations) # create array for cost value per iteration accuracy = np.zeros(iterations) # create array for achieved training accuracy per iteration # load and format datasets data = pd.read_csv('/kaggle/input/training-set/train.csv') # import *.csv using pandas data = np.array(data) # convert data into numpy array m, n = data.shape # m = number of samples, n= number of input neurons + label np.random.shuffle(data) # shuffle lines before splitting into test and training sets # create test dataset and reformat data_test = data[0:test_size].T # take 1st test_size rows and transpose Y_test = data_test[0] # 1st row contains labels X_test = data_test[1:n] / 255. # remaining rows contain normalized gray values # create training dataset and reformat data_train = data[test_size:m].T # treat training set in the same way Y_train = data_train[0] X_train = data_train[1:n] / 255. # define one hot representation of label vector def one_hot(Y): one_hot_Y = np.zeros((Y.max() + 1, Y.size)) # create zero matrix w/ shape (number of categories) x (number of samples) one_hot_Y[Y, np.arange(Y.size)] = 1 # set 1 in row Y for every column (from 0 to number of samples -1) return one_hot_Y Oh_Y_train = one_hot(Y_train) # calculate one_hot_y and set as a global variable Oh_Y_test = one_hot(Y_test) #define the neural network # define initial model parameter def init_params(): for l in range (1,L): params['W'+str(l)] = np.random.normal(size=(layer_dim[l], layer_dim[l-1])) * np.sqrt(1./layer_dim[l-1]) # create array of dim(layer_n+1) x dim(layer_n) with values around 0 and standard deviation of 1 scaled to sqrt(1/# neurons layer_n) params['b'+str(l)] = np.random.normal(size=(layer_dim[l], 1)) * np.sqrt(1./layer_dim[l]) # create array of dim(layer_n+1) x 1 with values around 0 and standard deviation of 1 scaled to sqrt(1/# neurons layer_n+1) return params # define non linear activation functions # define Rectified Linear Unit function for activation of inner layers def ReLU(Z): return np.maximum(Z, 0) # define softmax function to calculate propability for last layer (logistic regression) def softmax(Z): A = np.exp(Z) / sum(np.exp(Z)) return A # define forward propagation def forward_prop(X, params): cache['A0'] = X for l in range (1,LL): # calculate inner layers cache['Z'+str(l)] = params['W'+str(l)].dot(cache['A'+str(l-1)]) + params['b'+str(l)] # Z = W x A + b cache['A'+str(l)] = ReLU(cache['Z'+str(l)]) # A = ReLu(Z) cache['Z'+str(LL)] = params['W'+str(LL)].dot(cache['A'+str(LL-1)]) + params['b'+str(LL)] # calculate outer layer cache['A'+str(LL)] = softmax(cache['Z'+str(LL)]) # A_Last_Layer = softmax (Z_Last_Layer) return cache['A'+str(LL)], cache # define derivative of Rectified Linear Unit function def ReLU_deriv(Z): return Z > 0 # return all Z values > 0 , else return 0 # define backward propagation def backward_prop(params, cache, Oh_Y): grads['Z'+str(LL)] = cache['A'+str(LL)] - Oh_Y # calculate change of Z_LastLayer, Note: repeat {Z = Z -alpha * dL/dZ} for gradient descent, and dL/dA * dA/dZ = A-Y when using L = -(Y * ln(A) + (1-Y) * ln(1-A)) as a cost function, and softmax to activate A grads['W'+str(LL)] = 1 / m * grads['Z'+str(LL)].dot(cache['A'+str(LL-1)].T) # calculate change of W_LL , dW_l = 1/m * dZ_l * A_(l-1).T grads['b'+str(LL)] = 1 / m * np.sum(grads['Z'+str(LL)]) # calculate change of b_LL, dB_l = 1/m * sum(dZ_l) for l in range (LL-1,0,-1): # calculate change of inner layers grads['Z'+str(l)] = params['W'+str(l+1)].T.dot(grads['Z'+str(l+1)]) * ReLU_deriv(cache['Z'+str(l)]) # calculate change of Z_l = W_(l+1).T x dZ_(l+1) * dReLu(Z_l) grads['W'+str(l)] = 1 / m * grads['Z'+str(l)].dot(cache['A'+str(l-1)].T) # calculate change of W_l grads['b'+str(l)] = 1 / m * np.sum(grads['Z'+str(l)]) # calculate change of b_l return grads # define update of parameters def update_params(params, grads, alpha): for l in range (1,L): params['W'+str(l)] = params['W'+str(l)] - alpha * grads['W'+str(l)] # gradient descent W_l, W_l = W_l - alpha * dW_l params['b'+str(l)] = params['b'+str(l)] - alpha * grads['b'+str(l)] # gradient descent b_l, b_l = b_l - alpha * dB_l return params # define evaluation of model # define prediction function def get_predictions(AL): return np.argmax(AL, 0) # position of biggest value in Activated_LastLayer is predicted number # define calculation of accuracy def get_accuracy(predictions, Oh_Y): return np.sum(predictions == np.argmax(Oh_Y, 0)) / Oh_Y.shape[1] # sum up matches and normalize to samplesize # define calculation of cost function def compute_cost(AL, Oh_Y): cost = (-1 / m) * np.sum(np.multiply(Oh_Y, np.log(AL))) return (cost) # define cost and accuracy plot def plot_graph(x_values, cost, accuracy): fig, axs = plt.subplots(1, 2, figsize=(9, 3)) axs[0].set_xlabel('iteration') axs[0].set_ylabel('cost [a.u.]') axs[0].plot(x_values,cost) axs[1].set_xlabel('iteration') axs[1].set_ylabel('accuracy [%]') axs[1].set_ylim(0,100) axs[1].plot(x_values,accuracy) plt.show()
Nachdem alle notwendigen Funktionen definiert sind, starten wir das Training.
#start the training start = time.time() params = init_params() # initialize W and b for i in range(iterations): AL, cache = forward_prop(X_train, params) # perform forward propagation cost[i] = compute_cost(AL,Oh_Y_train) # calculate cost accuracy[i] = get_accuracy(get_predictions(AL), Oh_Y_train) # get average accuracy for all samples grads = backward_prop(params, cache, Oh_Y_train) # perform backward propagation params = update_params(params, grads, training_rate) # update model # print out cost and accurracy if i % 20 == 0: print("Iteration: "+str(i)) print("Cost :"+str(cost[i])+" Accuracy: "+str(accuracy[i])) end = time.time() duration = end - start # training duration in seconds plot_graph(list(range(0,iterations)), cost, accuracy*100) # plot cost and accuracy vs iterations
Bereits nach 200 Iterationen sättigt die Genauigkeit merklich. Wir erreichen eine Genauigkeit von ~85%.
Wir sichern das Model,
# Save model now = datetime.now() date = now.strftime('%Y-%m-%dT%H-%M-%S') + ('-%02d' % (now.microsecond / 10000)) name = "" for i in range (len(layer_dim)): name += str(layer_dim[i])+"_" name = date+"_dim_"+name+"samples_"+str(m)+"_rate_"+str(training_rate).replace('.', '')+"_iterations_"+str(iterations) np.save(name, params) # save to ./kaggle/working
und können es wieder laden.
# Load model model = np.load(name +".npy",allow_pickle='TRUE').item() # define how to predict label for a certain sample def make_prediction(X, model): AL ,_ = forward_prop(X, model) prediction = get_predictions(AL) return prediction # define how to test prediction vs label and image def test_prediction(index, model): current_image = X_test[:, index, None].reshape((28, 28)) * 255 prediction = make_prediction(X_test[:, index, None], model) label = Y_test[index] print("Prediction: ", prediction) print("Label: ", label) plt.gray() plt.imshow(current_image, interpolation='nearest') plt.show() return
Dann berechnen wir die für den Testdatensatz erzielte Genauigkeit. Alternativ könnten wir uns auch einige Testbilder und die aus dem Model bestimmte Kategorie anzeigen lassen.
# calculate overall accuracy of model for test dataset test_accurracy = get_accuracy(make_prediction(X_test, model), Oh_Y_test) print(test_accurracy)
Für den obigen Durchlauf erzielen wir auch für den nicht trainierten Testdatensatz eine Genauigkeit von fast 89%. Um die Performance der untersuchten Konfiguration zu dokumentieren, können wir abschließend noch ein log file speichern.
# write log file for training session log = ["date: "+date+"\n", "training samples: "+str(m)+"\n", "training duration: "+str(round(duration/60, 1))+" minutes\n", "training accuracy: "+str(round(accuracy[iterations-1]*100, 1))+" %\n", "test samples: "+str(test_size)+"\n", "test accurracy: "+str(round(test_accurracy*100, 1))+" %\n"] with open("log_"+str(name)+".txt", 'w') as f: f.writelines(log)
Viel Spaß beim selber programmieren...
Keine Kommentare:
Kommentar veröffentlichen