Eine selbstgebaute Baseball Wurfmaschine

Unser Sohn ist ein großer Baseball Fan. In den letzten Monaten sind daher ein Batting Cage fürs Schlagtraining im Carport entstanden,  sowie ein 'Screen' um die Bälle sicher zuwerfen zu können, und nun auch eine kleine Pitching Maschine.

Für die Wurfmaschine habe ich einen kleinen 230 V Elektromotor auf ein Holzbrettchen geklemmt. Der Motor treibt eine Kunststoffrolle mit einem Durchmesser von 10 cm, und läuft mit 6000 rpm, so dass sich theoretisch Ballgeschwindigkeiten von 113 km/h oder ca. 70 mph erzielen lassen. Dies entspricht zwar nicht den Geschwindigkeiten mit denen Profis der MLB werfen, ist aber ziemlich genau die Wurfgeschwindigkeit eines trainierten 14-Jährigen. Um den Anpressdruck der Rolle auf den Ball einzustellen, habe ich Kunststoffplatten unter den Motor geklemmt, bis die Bälle gut ausgeworfen werden.

Ein Abflussrohr mit 75 mm Durchmesser sorgt dafür das die Bälle in die gewünschte Richtung geschossen werden. Wir verwenden weiße Kunststoff Wiffle Bälle von Franklin. 

Diese passen genau in das Rohr und sind leicht deformierbar. Damit die Bälle am Aufsetzpunkt der Rolle nicht verkanten, habe ich das Abflussrohr seitlich geschlitzt. Die Maschine ist auf ein kleines Tischchen hinter den 'Screen' geklemmt.

Und so sieht die Maschine dann in Aktion aus.

Viel Spaß beim selber Basteln...

Erzeugung und Verteilung von Quantenschlüsseln

Am 16.August 2024 hat eine Falcon 9 Rakete der Firma Space X den Kleinsatelliten QUBE in den Orbit gebracht.

Mit QUBE können sogenannte Quantenschlüssel zwischen Satelliten und Bodenstationen ausgetauscht werden, um damit abhörsichere Kommunikation zu realisieren. 

Durch den Einsatz mehrerer solcher Kleinsatelliten könnte eine weltweite, sichere Kommunikation möglich werden, ohne auf große Glasfasernetze zurückgreifen zu müssen.  

Der QUBE-Satellit verwendet das BB84-Protokoll für die Quantenschlüsselverteilung oder Quantum Key Distribution (QKD). Dieses Quantenschlüsselverteilungsschema wurde bereits 1984 von Charles Bennett und Gilles Brassard entwickelt, und basiert auf der Polarisation von Photonen und der Tatsache, dass jede Messung eines Quantenzustands diesen verändert.

Der erste Schritt in BB84 ist eine Quantenübertragung. Der Sender Alice erstellt ein zufälliges Bit (0 oder 1), und wählt dann zufällig eine ihrer beiden Basen (in diesem Fall geradlinig oder diagonal) aus, um es zu übertragen. 

Anschließend bereitet sie einen Photonenpolarisationszustand vor, der sowohl vom Bitwert als auch von der Basis abhängt, wie in der obenstehenden Tabelle gezeigt. So ist z.B. eine 0 in der geradlinigen Basis (+) als vertikaler 0°-Polarisationszustand kodiert, und eine 1 ist in der diagonalen Basis (x) als 135°-Zustand. 

Alice sendet dann ein einzelnes Photon, wobei Alice den Zustand, die Basis und die Zeit jedes gesendeten Photons aufzeichnet.

Da der Empfänger Bob die Basis nicht kennt, in der die Photonen kodiert wurden, kann er nur zufällig eine Basis auswählen, in der er messen möchte, entweder geradlinig oder diagonal. Er tut dies für jedes Photon, das er empfängt, und zeichnet ebenfalls die Zeit, die verwendete Messgrundlage und das Messergebnis auf. 

Spielerisch lassen sich die Vorgänge gut mit Quantum Flytrap, einer No-Code-IDE für Quantencomputing, überprüfen.

Der Screenshot oben zeigt ein Virtual Lab, in dem Alice und Bob ihre Basiszustände mit Faraday Rotatoren einstellen. In der Simulation lässt sich, unter Anderem, die Zahl der durchzuführenden Messungen einstellen, und die Ergebnisse können als *.csv Datei exportiert werden.

Nachdem Bob alle Photonen gemessen hat, kommuniziert er mit Alice über einen öffentlichen klassischen Kanal. Alice sendet die Basis, nicht aber die Bits, in der jedes Photon eingesendet wurde, und Bob die Basis, in der jedes Photon gemessen wurde. Beide verwerfen Photonenmessungen (Bits), bei denen Bob eine andere Basis als Alice verwendet hat. Diese sind im Beispiel oben rot markiert. 

Um zu überprüfen, ob ein Lauscher vorhanden ist, vergleichen Alice und Bob nun eine vorgegebene Teilmenge ihrer verbleibenden Bitzeichenfolgen. Wenn eine dritte Partei (normalerweise als Eve bezeichnet, für "Lauscherin") Informationen über die Polarisation der Photonen erhalten hat, führt dies zu Fehlern in Bobs Messungen, und die gemessenen Bits stimmen nicht mit den gesendeten überein.

War kein Lauscher vorhanden bilden die nicht verglichenen Bits den Schlüssel. 

Die eigentliche Nachricht kann anschließend z.B. mit Hilfe der One Time Pad Verschlüsselung (OTP) codiert werden. Dafür ist es notwendig, dass ein Schlüssel verwendet wird, der mindestens so lang wie die Nachricht ist.

Nachricht und Schlüssel werden miteinander durch eine Exklusiv Oder Operation verknüpft. Das OTP ist informationstheoretisch sicher und kann nachweislich nicht gebrochen werden.

Viel Spaß beim Spielen mit der Quantenkryptographie...

Ein Untergestell für die Drechselbank

Meine Frau hat sich zum Geburtstag eine Drechselbank gewünscht. Zum Einstieg in dieses neue Hobby fiel meine Wahl auf die PDB 100 A1 Drechselbank von Parkside. Die Maschine verfügt über eine Motorleistung von 400W, und durch Umlegen eines Keilriemens lassen sich Drehzahlen zwischen 890, 1260, 1760 und 2600 U/min wählen. Das Bankbett lässt sich von 50 cm auf 1 m verlängern, so dass auch größere Werkstücke verarbeitet werden können.

Ich habe die Maschine auf eine Arbeitsplatte aus Eiche montiert, und dieses auf ein kleines Untergestell geschraubt, das ich an der Wand befestigt habe. An der rechten Seitenwand ist eine Halterung für die Drechseleisen vorgesehen. 

Im Zuge des weiteren Ausbaus habe ich noch ein horizontales Brett in das Untergestell eingefügt, um eine zusätzliche Ablage für Werkzeug zu haben.

Das Bankbett der Maschine habe ich durch 3 mm starkem Flacheisen verstärkt. Das Flacheisen ist an mehreren Stellen mit der Maschine verschraubt, und verhindert ein Durchbiegen des gefalzten Bleches beim festen Einspannen der Werkstückes mit dem Reitstock. 

Zum Schärfen der Drechseleisen habe ich eine Nassschleifmaschine TIGER 2000S von Scheppach erworben.

Die Maschine wurde mit etlichem Zubehör geliefert. Um aber insbesondere die für das Langholzdrechseln wichtige Schruppröhre schärfen zu können, habe ich zusätzlich in die Drechselmesser-Schleifführung TWSTGJ von Triton investiert. Diese erlaubt durch Kippen der Rundröhren um einen fixen Auflagepunkt den sogenannten Fingernagelschliff.

Viel Spaß beim selber Drechseln...

Eine sprachgesteuerte Balleinwurfmaschine für den Kickerkasten

Ein Tischfußball begeisterter Freund hat mich auf die Idee gebracht eine Einwurfmaschine für den Kickerkasten zu bauen. Der Ball soll auf einen Sprachbefehl hin eingeworfen werden, und die Maschine soll an verschiedenen Stellen im Spielfeld platziert werden können.

Im Bild oben ist meine Realisierung zu sehen. Eine gedruckte Scheibe, mit einer Tasche für den Ball, dreht sich in einem Kunststoff Topf. Der Ball wird während einer Umdrehung der Scheibe aufgenommen und fällt dann dann durch ein Loch im Boden des Topfes. Eine röhrenförmige Barriere verhindert, dass mehr als ein Ball durch das Loch fallen kann. In der Ansicht von unten erkennt man, dass der Ball durch ein gewinkeltes Stück Abflussrohr mit dem Durchmesser von 50 mm nach vorne ausgeworfen wird.

Als Antrieb habe ich einen 12V Schneckengetriebemotor mit 30rpm verbaut.

Laut Herstellerangaben verfügt dieser über ein maximales Drehmoment von 10 kg x cm, und dreht sich daher auch in einem mit Bällen voll gefüllten Topf gut.

Die gedruckte Scheibe mit der Tasche besteht aus zwei Teilen. Der untere Teil

enthält die Aufnahme für die Motorachse und Abstandshalter. Der obere Teil besteht aus einem planen Deckel, der sich auf die Abstandshalter schrauben lässt. 

Der Kunststofftopf mit den Bällen ist auf laminierte Sperrholzplatten der Stärke 18mm geschraubt, die sich auf die Stangen des Kickertisches stellen lassen. 

Im Bild oben ist rechts unten bereits die Steuerungseinheit der Maschine und ein Konferenz Mikrofon zu erkennen. Das gedruckte Gehäuse der Steuerung besteht aus Deckel und Boden,

und enthält neben einem Raspberry Pi 3 Modell B+ auch ein Relais, das mit 3V ansteuerbar ist, und Ströme bis zu 10A schalten kann. 

Das im Bild oben erkennbare USB Dongle Mikrofon mit einer Empfindlichkeit von -58 dB laut Herstellerangabe, habe ich später durch ein Konferenz-Mikrofon mit einer Empfindlichkeit von -45 dB ersetzt, um die Spracheingabequalität zu verbessern.

Für die Sprachsteuerung habe ich das Open-Source Software Toolkit VOSK von Aphacephei genutzt. Es unterstützt mehr als 20 Sprachen und Dialekte, darunter Englisch und Deutsch. Die kleinen Sprachmodelle, die sich auf dem Github repository von Alphacephei finden, laufen offline auch auf Geräten wie dem Raspberry Pi.  

Eine Einführung in des  Toolkit findet sich auf der Alphacephei website. Nach der Installation mit

pip3 install vosk

habe ich das Sprachmodell 'vosk-model-small-de-0.15.zip' durch

git clone https://github.com/alphacep/vosk-api

cd vosk-api/python/example

wget https://alphacephei.com/vosk/models/vosk-model-small-de-0.15.zip

unzip vosk-model-small-de-0.15.zip

mv vosk-model-small-de-0.15/ model

aus dem repository geladen und entpackt.

Die Ablaufsteuerung erfolgt in Python. Der folgende Code wurde mit dem Syntax Highlighter highlight.hohli von Anton Shevchuk für den Blog aufbereitet, und mit 

<div style="background: rgb(255, 255, 255); overflow: auto; width: auto;"></div> 

ergänzt, um horizontales scrolling des Abschnitts zu gewährleisten.

import queue
import json
import sounddevice as sd
import vosk
import time
import gpiozero
import sys

q = queue.Queue()
STARTCODEs = ["einwurf", "ein wurf", "ein uhr", "ein ver", "nun eingabe", "eingabe", "ein einwurf", "ein", "ein waffen", "nun ein nur", "nun einwurf", "eingabe", "mehr", "auch ein mensch", "einwirken", "nun ein wolf", "eingriff"]

def callback(indata, frames, time, status):
    #"""This is called (from a separate thread) for each audio block."""
    if status:
        print(status, file=sys.stderr)
        pass
    q.put(bytes(indata))

if __name__ == '__main__':

    #Relais initialisieren
    red_led = gpiozero.LED(4)
    green_led = gpiozero.LED(17)

    # Speech Objekt erstellen
    model = vosk.Model('/home/pi/vosk-api/python/example/model')
    
    #Stream vom Mikrofon öffnen
    with sd.RawInputStream(samplerate=44100, blocksize=16000, device=None,dtype='int16',
                            channels=1, callback=callback):
        #Sternchen zeichnen
        #print('*' * 80)
        #Programmstart mit roter LED anzeigen
        red_led.on()
        
        # Aktivierung der vosk Spracherkennung mit Übergabe des geladenen Models. Übersetze das Gesprochene in Text.
        rec = vosk.KaldiRecognizer(model, 44100)
        while True:
            # Daten aus der Queue ziehen
            data = q.get()
            #print("Bitte sprechen!")
            if rec.AcceptWaveform(data):
                # wandle das gesprochene Wort in jason (JavaScript Object Notation) 
                res = json.loads(rec.Result())
                # und gib es aus
                print(res)
                # wenn der Aktivierungscode herausgehört wurde, wird das Relais geschaltet
                for wort in STARTCODEs:
                    if wort == res['text']:
                        green_led.on()
                        red_led.off()
                        # Relais 1 Sekunden anziehen
                        time.sleep(1)
                        green_led.off()
                        red_led.on()

            else:
                pass

Nach Laden der vosk Bibliothek wird im Hauptteil des Programms beim Start die rote Status LED gesetzt. Mit der 'sounddevice' Funktion 'RawInputStream' werden Aufnahmestückchen einer bestimmten Länge in eine 'queue' geschrieben.  Werden die Stückchen vom Modell als 'waveform' akzeptiert, werden sie transkribiert und mit den in 'STARTCODES' definierten Begriffen verglichen. Wenn der richtige Aktivierungscode herausgehört wurde, zieht das Relais für 1.3 Sekunden an, und die grüne LED wird gesetzt. Die Drehscheibe vollführt etwas mehr als eine Umdrehung, und wirft einen Ball aus. Danach erlischt die grüne LED, und die rote LED wird wieder gesetzt. Das System ist bereit für  das nächste Aktivierungswort. 

Damit das Programm nach dem Booten direkt startet, muss mit 'systemd' ein Service gestartet werden. Eine detaillierte Anleitung wie ein solcher 'Service' erzeugt werden kann findet sich z.B. hier.

Mein '/etc/systemd/system/AutoRunEinwurf_Service' file sieht so aus.

Zum bequemen Programmieren via Remote Desktop vom Laptop aus, habe ich diesesmal VNC statt X11 benutzt. Der VNC Server wird von Raspberry Pi OS standardmäßig unterstützt, muss aber aktiviert werden. Um des zu tun kann man raspi-config auf der Befehlszeile mit

sudo raspi-config

öffnen. Anschließend muss man zu 'Schnittstellenoptionen' navigieren, VNC auswählen und die bestätigen.

Als Client auf dem Laptop habe ich tigervnc heruntergeladen. Nach der Installation kann man sich wie gewohnt durch Eingabe der IP Adresse des Raspberry

auf den Desktop des kleinen Rechners einloggen.

Zu guter Letzt möchte ich die Maschine noch in Aktion zeigen.

Viel Spaß beim selber Basteln...

Die Maxwellgleichungen anschaulich verstehen

Die Maxwellgleichungen beschreiben die grundlegenden physikalischen Gesetze im Zusammenhang mit Elektrizität und Magnetismus. 

James Clerk Maxwell

Die mathematische Formulierung  erscheint zunächst schwierig, die Gleichungen lassen sich aber anschaulich verstehen.

\[\begin{flalign} & (1) \qquad \operatorname{div}\mathbf{E} = \frac{\rho}{\varepsilon_0} \\ & (2) \qquad \operatorname{div} \mathbf{B} = 0 \\ & (3) \qquad \operatorname{rot} \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \\ & (4) \qquad \operatorname{rot} \mathbf{B} = \mu_0 \left(\mathbf{J} + \varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \end{flalign}\]

Bertrachten wir zunächst die erste Gleichung, das sogenannte Gaußsche Gesetz.

  \[\operatorname{div} \mathbf{E} = \frac{\rho}{\varepsilon_0}\]

Die Divergenz \(\operatorname{div}\) eines Vektorfeldes \(\mathbf{E}\) ist dabei als Skalarprodukt des Nabla Operators \( \nabla = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z}\right)\) mit dem Feld \(\mathbf{E}\), also als \(\nabla \cdot \mathbf{E}\) definiert.

Entfernen wir uns im zwei dimensionalen Raum einen kleinen Schritt von der Position \(\left(x_0,y_0\right)\), so zeigt das Vektorfeld in der Regel in eine andere Richtung. Der rote Vektor im Bild oben ist der Differenzvektor.

Betrachten wir alle Differenzvektoren in der Umgebung von \(\left(x_0,y_0\right)\), und mitteln über die Skalarprodukte aus Abstandsvektor und Differenzvektor, so ergibt sich im Grenzübergang zu sehr kleinen Abständen \(\left( \partial / \partial x, \partial / \partial y\right) \cdot \left(E_x, E_y\right) \). Da sich das Skalarprodukt aus dem Produkt der Beträge der Vektoren mal dem Cosinus des Zwischenwinkels errechnet, wird es maximal wenn Abstandsvektor und Differenzvektor in die gleiche Richtung zeigen. Der Ausdruck \(\nabla\cdot \mathbf{E}\) beschreibt also wieviel 'Feld' aus dem Punkt \(\left(x_0,y_0\right)\) gemittelt über alle Richtungen 'herausfließt'. Das Gaußsche Gesetzt sagt nun, dass die 'Menge des herausfließenden elektrischen Feldes' proportional zur Punktladungsdichte ist. Punktladungen sind also die Quellen und Senken des elektrischen Feldes.

Die zweite Gleichung 

\[\operatorname{div} \mathbf{B} = 0\]

wird häufig auch als Gaußsches Gesetz des Magnetismus bezeichnet. Die Divergenz der magnetischen Flussdichte \(\mathbf{B}\) ist Null. Es gibt keine Quellen oder Senken, d.h. keine magnetischen Monopole, die analog zu elektrostatischen Ladungen wirken würden. Das magnetische Feld verhält sich wie eine nicht komprimierbare Flüssigkeit. 

Die dritte Gleichung

   \[\operatorname{rot} \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}\]

wird als Faradaysches Induktionsgesetzt bezeichnet.

Als Rotation \(\operatorname{rot}\) des Vektorfeldes \(\mathbf{E}\) bezeichnet man dabei das Vektorprodukt des Nabla Operators \( \nabla = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z}\right)\) mit dem Feld \(\mathbf{E}\), also \(\nabla \times \mathbf{E}\).

Im Gegensatz zum Skalarprodukt, ist das Vektorprodukt ein Maß dafür, ob zwei Vektoren zu einander senkrecht stehen. Sein Betrag errechnet sich als Produkt der Beträge, mal dem Sinus des Zwischenwinkels. Der Betrag des Vektorprodukts wird also für einen Zwischenwinkel von 90° maximal. Betrachten wir wieder die Differenzvektoren in einigem Abstand von \(\left(x_0,y_0\right)\), so entspricht die Rotation \(\operatorname{rot}\) dem Grenzwert des Mittelwertes aus Abstandsvektoren \(\times\) Differenzvektoren für kleine Abstände im Übergang zum Punkt \(\left(x_0,y_0\right)\). Das Induktionsgesetz besagt nun, dass eine zeitliche Änderung des Magnetfeldes bzw. der magnetischen Flussdichte zu einer Rotation des elektrischen Feldes führt. Anders ausgedrückt resultiert also die zeitliche Änderung des durch eine Leiterschleife eingeschlossenen magnetischen Flusses in einer Spannung an den Enden der Leiterschleife. Die induzierte Spannung wirkt der Änderung entgegen. Das Bild unten von Michael Lenz veranschaulicht dies.

Die vierte Gleichung

\[ \operatorname{rot} \mathbf{B} = \mu_0 \left(\mathbf{J} + \varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \]

wird als Ampèresches Gesetz mit Maxwells Ergänzung bezeichnet. Der erste Term auf der rechten Seite der Gleichung besagt, dass die Rotation des Magnetfeldes bzw. die magnetische Flußdichte \( \mathbf{B}\) proportional zum elektrischen Stromfluss bzw. der Stromdichte \( \mathbf{J}\) ist. Ist das elektrische Feld \(\mathbf{E}\) zeitlich konstant, können wir den zweiten Term vernachlässigen. Für einen von Gleichstrom durchflossenen Leiter gilt also das untenstehende Bild.

Ändert sich das elektrische Feld \(\mathbf{E}\), ist die Rotation zusätzlich proportional zur zeitlichen Änderung des elektrischen Feldes. Diese Änderung wird häufig auch Verschiebungsstrom genannt, und stellt die sogenannte Maxwell Ergänzung dar.

Aus der Maxwell Ergänzung im Wechselstromfall

\[ \operatorname{rot} \mathbf{B} = \mu_0 \left(\varepsilon_0\frac{\partial \mathbf{E}}{\partial t} \right) \]

lässt sich durch erneutes Anwenden des Rotationoperators recht einfach die Wellengleichung

$$\Delta \mathbf{E} = \mu_0 \varepsilon_0 \frac{\partial^2 \mathbf{E}}{\partial t^2} $$

herleiten, \(\Delta\) bezeichnet dabei den sogenannten Laplace Operator. Eine Lösung dieser Differentialgleichung sind ebene Wellen der Form

$$\mathbf{E} = E_0 f(k \cdot \mathbf{x} - ct)$$

Sie beschreiben die Ausbreitung einer elektromagnetischen Welle im Raum.

Viel Spaß beim selber Rechnen und Verstehen...

Ein neuronales Netzwerk von Grund auf programmieren

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...