Die Jupyter-Notebooks zur Lehrveranstaltung finden Sie im zugehörigen GitHub-Repository.

Vektoren, Matrizen und Vektorisierung

3. Vektoren, Matrizen und Vektorisierung in Python¶

In dieser Einheit werden wir uns einerseits einige Werkzeuge für das Rechnen mit Vektoren und Matrizen in Python ansehen. Andererseits möchte ich Ihnen das Prinzip der Vektorisierung klar machen (und nahelegen).

Hier zunächst die Imports für heute:

In [24]:

%matplotlib inline
import matplotlib.pyplot as plt # für plotting, wie gewohnt
from matplotlib import image    # zum Einlesen und arbeiten mit Bilddateien

import numpy as np              # für numerische Aktionen mit Arrays, wie gewohnt

from tqdm import tqdm           # zum Erzeugen eines Fortschrittsbalkens

import os                       # Funktionen zum OS

3.1 Vektoren¶

Vektoren sind uns allen soweit bekannt (zumindest sollen sie das sein). In diesem Abschnitt geht es daher in erster Linie darum, wie Sie bestimmte gewohnte Dinge mit Vektoren tun können. Wir verwenden hier grundsätzlich die Package NumPy.

In [25]:

# Ein Vektor ist in Python im Prinzip ein 1D-Array
a_vector = np.array([2, 3, -1])

# Zwei Vektoren lassen sich recht einfach addieren und subtrahieren
another_vector = np.array([4, -2, 1])

print("Die Summe: ", a_vector + another_vector)

print("Die Differenz: ", a_vector - another_vector)

print("Das punktweise Produkt: ", a_vector * another_vector)

print("Der punktweise Quotient: ", a_vector / another_vector)

print("Punktweises Potenzieren: ", a_vector.astype(float) ** another_vector)

Die Summe:  [6 1 0]
Die Differenz:  [-2  5 -2]
Das punktweise Produkt:  [ 8 -6 -1]
Der punktweise Quotient:  [ 0.5 -1.5 -1. ]
Punktweises Potenzieren:  [16.          0.11111111 -1.        ]

Was ist aber mit dem Skalar-Produkt $\vec{a}\cdot\vec{b}$ oder dem Kreuz-Produkt $\vec{a}\times\vec{b}$? Eine Möglichkeit ist es, NumPy Funktionen dafür zu verwenden und die entsprechenden Ausdrücke zu basteln. Oder, man bedient sich fertiger Funktionen, die ebenfalls in NumPy zu finden sind.

In [26]:

# bastel-Variante für das Skalarprodukt
print(np.sum(a_vector * another_vector))

# built-in-Variante für das Skalarprodukt
print(np.dot(a_vector, another_vector))

# built-in-Variante für das Kreuzprodukt
print(np.cross(a_vector, another_vector))

1
1
[  1  -6 -16]

3.2 Matrizen und Lineare Algebra¶

Matrizen sind in Python im Prinzip 2-dimensionale Arrays. Man gibt einem NumPy-Array einen Zeilen- und einen Spalten-Index und kann damit dann verschiedene Dinge tun (das meiste davon geht natürlich auch schon in einer Dimension, ist aber in 2D anschaulicher), z.B.:

Die Matrix transponieren
Einen Teil der Matrix herausschneiden
Elemente an einer bestimmten Position verändern, hinzufügen, löschen
Matrizen miteinander (also auch mit Vektoren) multiplizieren
Eigenschaften der Matrix berechnen, wie z.B. die Determinante, den Rang, etc.
Operationen aus der Linearen Algebra ausführen, wie z.B. Eigenwerte bestimmen

Hier eine Auswahl der Dinge auf dieser Liste:

In [27]:

# definiere eine Matrix
a_matrix = np.array([[2, 3, 5, -3], [1, 3, -1, 4], [1, 2, -3, 2], [3, 4, -2, 1]])
a_matrix

Out[27]:

array([[ 2,  3,  5, -3],
       [ 1,  3, -1,  4],
       [ 1,  2, -3,  2],
       [ 3,  4, -2,  1]])

In [28]:

# verschiedene Slices/Teile
print("Erste Zeile: ", a_matrix[0, :])

print("Erste Spalte: ", a_matrix[:, 0])

print("Letzte Zeile: ", a_matrix[-1, :])

print("Letzte Spalte: ", a_matrix[:, -1])

print("Element links oben: ", a_matrix[0, 0])

print("Element an bestimmter Position: ", a_matrix[2, 3])

print("Elemente im Viereck innen: \n", a_matrix[1:3, 1:3])

print("Jedes 2. Element der 2. Zeile: ", a_matrix[1, ::2])

print("Dritte Spalte, von unten nach oben: ", a_matrix[::-1, 2])

Erste Zeile:  [ 2  3  5 -3]
Erste Spalte:  [2 1 1 3]
Letzte Zeile:  [ 3  4 -2  1]
Letzte Spalte:  [-3  4  2  1]
Element links oben:  2
Element an bestimmter Position:  2
Elemente im Viereck innen: 
 [[ 3 -1]
 [ 2 -3]]
Jedes 2. Element der 2. Zeile:  [ 1 -1]
Dritte Spalte, von unten nach oben:  [-2 -3 -1  5]

In [29]:

# transponierte Matrix
np.transpose(a_matrix)

Out[29]:

array([[ 2,  1,  1,  3],
       [ 3,  3,  2,  4],
       [ 5, -1, -3, -2],
       [-3,  4,  2,  1]])

In [30]:

# Elemente verändern
a_changed_matrix = a_matrix

a_changed_matrix[0, 1] = -4
second_row = a_changed_matrix[1]

print("Matrix geändert auf \n", a_changed_matrix)
print("Die zweite Zeile: ", second_row)

# ein Element an bestimmter Stelle in bestimmter Richtung herauslöschen
a_changed_matrix = np.delete(a_changed_matrix, 1, axis=0)

print("Zweite Zeile wurde gelöscht: \n", a_changed_matrix)

# Die gelöschte Zeile wieder einfügen
a_changed_matrix = np.insert(a_changed_matrix, 1, second_row, axis=0)

print("Zweite Zeile wieder eingefügt: \n", a_changed_matrix)

# eine Zeile anhängen
a_changed_matrix = np.append(a_changed_matrix, [second_row], axis=0)

print("Zweite Zeile nochmal angehängt: \n", a_changed_matrix)

# zwei Arrays kombinieren
a_changed_matrix = np.concatenate((a_changed_matrix, a_changed_matrix), axis=1) # versuche auch axis=0

print("Die veränderte Matrix, kombiniert mit sich selbst: \n", a_changed_matrix)

# und, wichtig: Die Dimensionen der Matrix verändern (also das Array verformen)
a_changed_matrix = np.reshape(a_changed_matrix, (2, -1))

print("zwei Zeilen und wie viele Spalten? \n", a_changed_matrix)

Matrix geändert auf 
 [[ 2 -4  5 -3]
 [ 1  3 -1  4]
 [ 1  2 -3  2]
 [ 3  4 -2  1]]
Die zweite Zeile:  [ 1  3 -1  4]
Zweite Zeile wurde gelöscht: 
 [[ 2 -4  5 -3]
 [ 1  2 -3  2]
 [ 3  4 -2  1]]
Zweite Zeile wieder eingefügt: 
 [[ 2 -4  5 -3]
 [ 1  3 -1  4]
 [ 1  2 -3  2]
 [ 3  4 -2  1]]
Zweite Zeile nochmal angehängt: 
 [[ 2 -4  5 -3]
 [ 1  3 -1  4]
 [ 1  2 -3  2]
 [ 3  4 -2  1]
 [ 1  3 -1  4]]
Die veränderte Matrix, kombiniert mit sich selbst: 
 [[ 2 -4  5 -3  2 -4  5 -3]
 [ 1  3 -1  4  1  3 -1  4]
 [ 1  2 -3  2  1  2 -3  2]
 [ 3  4 -2  1  3  4 -2  1]
 [ 1  3 -1  4  1  3 -1  4]]
zwei Zeilen und wie viele Spalten? 
 [[ 2 -4  5 -3  2 -4  5 -3  1  3 -1  4  1  3 -1  4  1  2 -3  2]
 [ 1  2 -3  2  3  4 -2  1  3  4 -2  1  1  3 -1  4  1  3 -1  4]]

In [31]:

# zwei Matrizen miteinander multiplizieren
a_matrix_squared = np.matmul(a_matrix, a_matrix)
a_matrix_squared

Out[31]:

array([[ -4, -22,   5, -15],
       [ 16,  19,  -3,  11],
       [  7,   4,   8,   1],
       [ 11,   0,  15,   4]])

In [32]:

# Ein paar Dinge aus der Linearen Algebra

print("Rang der Matrix ", np.linalg.matrix_rank(a_matrix))

print("Determinante der Matrix ", np.linalg.det(a_matrix))

print("Eigensystem der Matrix: ")
print("Eigenwerte der Matrix: ", np.linalg.eig(a_matrix)[0])
print("Eigenvektoren der Matrix: \n", np.linalg.eig(a_matrix)[1])

Rang der Matrix  4
Determinante der Matrix  126.0
Eigensystem der Matrix: 
Eigenwerte der Matrix:  [ 3.82292406+2.38876617j  3.82292406-2.38876617j -2.32292406+0.89695077j
 -2.32292406-0.89695077j]
Eigenvektoren der Matrix: 
 [[-0.71322195+0.j         -0.71322195-0.j          0.10244504-0.22757172j
   0.10244504+0.22757172j]
 [ 0.40254467+0.3848397j   0.40254467-0.3848397j   0.47509046+0.19209634j
   0.47509046-0.19209634j]
 [ 0.09504103+0.18690531j  0.09504103-0.18690531j -0.10378635+0.36880975j
  -0.10378635-0.36880975j]
 [ 0.05505865+0.36629606j  0.05505865-0.36629606j -0.7268507 +0.j
  -0.7268507 -0.j        ]]

3.3 Beispiel: Lineare Regression¶

Um auch etwas konkret Sinnvolles mit Matrizen zu machen, greifen wir kurz auf die Optimierungsmethode der linearen Regression vor. Keine Sorge, wir sehen uns das später noch im Detail an. Einstweilen hier einmal nur der letzte Schritt der Anwendung.

Zunächst aber die wichtigsten Punkte: Die Lineare Regression leistet folgendes:

Man hat einen Datensatz, bei dem es eine bis mehrere (insgesamt $N$) unabhängige Variablen $x_i$ und eine abhängige Variable $y$ gibt.
Man möchte die Daten beschreiben, und zwar mit einer Funktion $f(x_i)$.
Bei der Linearen Regression setzt man diese Funktion linear an, also $f(x_i)=\sum_{i=0}^N a_i x_i$
Nimmt man dabei $x_0=1$ an, dann hat man den konstanten Term $a_0$ im Fit mit dabei.
Nun hat man einen Datensatz von $M$ Datenvektoren $\vec{x}^{(j)}$ zur Verfügung, also $j=1,\ldots,M$
Die Koeffizienten $a_i$ bestimmt man so, dass die Summe der Quadrate der Differenzen zwischen den $y^{(j)}$ und den $f(x_i)^{(j)}$ minimal ist (least-mean-squares fit).

Das sieht auf den ersten Blick verwirrend aus, aber worum es sich hier handelt, ist im Allgemeinen ein überbestimmtes System von linearen Gleichungen. Wir werden es uns außerdem anschaulich und leicht machen und nehmen eine einzige unabhängige Variable $x_1=x$ und die Konstante über $x_0=1$ mit.

In [33]:

# generieren wir erst einmal Daten
x_values = np.linspace(0, 10, 200)
fake_data = 0.5 * x_values + np.random.normal(size=len(x_values))

fig = plt.figure()

plt.scatter(x_values, fake_data, c="r", marker="x")

plt.show()

Vektoren, Matrizen und Vektorisierung in Python 2

Hierfür sollen wir nun den besten linearen Fit finden. Dafür gibt es zwar auch iterative Algorithmen wie gradient descent, aber wir werden hier eine geschlossene Formel für die Lineare Regression verwenden.

Dabei schreibt man den Parameter-Vektor der $a_i$ als $\vec{a}$, den Vektor der Funktionswerte $y^{(j)}$ als $\vec{y}$ und schließlich die $x$ Werte als Matrix $X$. Warum als Matrix? Das kommt daher, dass es einerseits $N$ unabhängige Variablen und andererseits $M$ Datenpunkte gibt, d.h. $X=x_i^{(j)}$. Konkret ist dabei in $X$ das $i$ der Zeilen- und $j$ der Spalten-Index.

Dann lautet die Lösung für die Koeffizienten folgendermaßen: $$\vec{a}=(X^TX)^{-1}\,X^T \vec{y}$$ Und das setzen wir nun im Code um, und plotten dann die entsprechende Gerade zu den Daten dazu:

In [34]:

# Konstuiere y und X
y_vec = fake_data

# für die Matrix X werden x-Werte und Einsen "nebeneinander" geschrieben
x_mat = np.hstack((x_values.reshape((len(x_values),1)), np.ones((len(x_values),1))))

# berechne a nach der Formel. Das "@" ist eine Kurzschreibweise für Matrixmultiplikation
a_vec = np.matmul(np.linalg.pinv(np.transpose(x_mat) @ x_mat), np.transpose(x_mat) @ y_vec)

# die Lösung ist ein Vektor aus a_1 und a_0
print("Die fit-Koeffizienten: ", a_vec)

# Erzeuge die Fit-Gerade als Einsetzen der x-Werte in ein Polynom mit den Koeffizienten a_1 und a_0
reg_data = np.polyval(a_vec, x_values)

# neue Figur
fig = plt.figure()

# plotte wieder die Daten
plt.scatter(x_values, fake_data, c="r", marker="x")

# und dazu die Gerade
plt.plot(x_values, reg_data, "b", label="Linear fit")

# die Legende
plt.legend(loc=0)

plt.show()

Die fit-Koeffizienten:  [0.48495158 0.03093417]

Vektoren, Matrizen und Vektorisierung in Python 3

3.4 Höherdimensionale Arrays¶

Arrays können im Prinzip auch mehr als nur 2 Dimensionen haben. Man muss nicht unbedingt in die tiefen irgendwelcher Theorien eintauchen, um so etwas zu finden. In der Bildbearbeitung (und auch z.B. im maschinellen Lernen an Bildern) hat man mit einem 2-dimensionalen Array zu tun, das noch drei weitere “Kanäle” hat. Im Daten-Array für das Bild wird das als dritte Dimension angelegt, sodass ein Bild-Array im Allgemeinen die Dimensionen (Breite in Pixel)x(Höhe in Pixel)x3 hat.

Wir werden uns hier ein Bild hernehmen, um ein paar Tests und einfache Manipulationen auszuprobieren. Sie werden sehen, dass vieles recht einfach ist und Ihnen bereits ein gutes Gefühl für den Umgang mit höherdimensionalen Daten gibt.

In [35]:

# lade ein Bild über matplotlib-Bibliothek "image"
# os.path.join zum OS-unabhängigen aneinanderfügen von Dateipfaden
the_picture = image.imread(os.path.join('data', 'pinky.jpg'))

# schauen wir uns das Bild an
plt.figure()

# der Befehl für die Darstellung eines Bildes
plt.imshow(the_picture)

plt.show()

In [36]:

# eigentlich ist das Bild aber "nur" ein Array, wie oben beschrieben
the_picture

Out[36]:

array([[[224, 220, 208],
        [224, 220, 208],
        [224, 220, 208],
        ...,
        [206, 211, 214],
        [205, 210, 213],
        [205, 210, 213]],

       [[226, 222, 210],
        [225, 221, 209],
        [225, 221, 209],
        ...,
        [207, 212, 215],
        [207, 212, 215],
        [207, 212, 215]],

       [[227, 223, 211],
        [227, 223, 211],
        [226, 222, 210],
        ...,
        [209, 214, 217],
        [209, 214, 217],
        [209, 214, 217]],

       ...,

       [[200, 189, 187],
        [200, 189, 187],
        [200, 189, 187],
        ...,
        [159, 150, 143],
        [157, 148, 141],
        [156, 147, 140]],

       [[197, 186, 184],
        [197, 186, 184],
        [197, 186, 184],
        ...,
        [161, 154, 146],
        [155, 148, 140],
        [148, 141, 133]],

       [[194, 183, 181],
        [194, 183, 181],
        [195, 184, 182],
        ...,
        [152, 145, 137],
        [154, 147, 139],
        [152, 145, 137]]], dtype=uint8)

In [37]:

# Die Dimension des Arrays:
np.shape(the_picture)

Out[37]:

(1512, 2016, 3)

In [38]:

# Wir können diese Werte auch so normieren, dass sie zwischen 0 und 1 liegen:
the_picture_normalized = the_picture/256


# jetzt können wir mit "matshow" statt "imshow" weitermachen
plt.matshow(the_picture_normalized)

Out[38]:

<matplotlib.image.AxesImage at 0x7fa960ee6be0>

Hier kommen nun einige einfache Manipulationen und Tests, die wir mit dem Bild machen können, indem wir einfach nur das Array verändern:

Spiegelung
Schwarz-Weiß
Nur einen Farbkanal darstellen
Aufhellen
Einen Filter über das Bild laufen lassen
- Unschärfe
- Beliebiger Filter

In [39]:

# gespiegelt, rechts-links

# Slicen aller Indices der Dimension 1 in umgekehrter Reihenfolge
plt.matshow(the_picture_normalized[:,::-1,:])

plt.show()

Vektoren, Matrizen und Vektorisierung in Python 6

In [40]:

# schwarz-weiß

# Hier brauchen wir eine Definition, wie man aus 
# RGB-Werten (also den drei Farbkanälen) einen
# Grayscale-Kanal macht. Die einfachste Möglichkeit dafür
# ist ein arithmetisches Mittel:
the_picture_normalized_bw = np.mean(the_picture_normalized, axis=2)

# danach hat das Array eine Dimension verloren und sieht seltsam aus
print(np.shape(the_picture_normalized_bw))

# deshalb müssen wir spezifisch eine Colormap (cmap) angeben
# die richtige für grayscale ist "gray" (siehe unten), aber es gibt
# viele teils lustige Möglichgkeiten: Greys, spring, magma, PRGn, seismic, hsv, Set2, prism
plt.matshow(the_picture_normalized_bw, cmap="gray") 

plt.show()

(1512, 2016)

In [41]:

# Nur ein Farbkanal

# Hier suchen wir uns also einen Kanal aus, z.B. R
the_picture_normalized_one_channel = the_picture_normalized[:,:,0]

# danach hat das Array wieder eine Dimension verloren
print(np.shape(the_picture_normalized_one_channel))

# deshalb müssen wir wieder eine Colormap (cmap) angeben
# die spezifischen hier sind "Reds", "Greens" und "Blues"
plt.matshow(the_picture_normalized_one_channel, cmap="Reds") 

plt.show()

(1512, 2016)

Vektoren, Matrizen und Vektorisierung in Python 8

In [42]:

# Alle drei Farbkanäle nebeneinander
plt.figure(figsize=(10, 5))

# Definiere alle nötigen Colormaps in einer Liste
colormaps = ["Reds", "Greens", "Blues"]

# Loop über den index der drei Farbkanäle
for ind_color in range(3):
    ax = plt.subplot(1, 3, ind_color+1)
    ax.matshow(the_picture_normalized[:,:,ind_color], cmap=colormaps[ind_color]) 

plt.show()

In [43]:

# Aufhellen
plt.figure(figsize=(10, 5))

# normales Bild zum Vergleich
ax = plt.subplot(1, 2, 1)
ax.matshow(the_picture_normalized) 

# aufgehelltes Bild daneben
ax = plt.subplot(1, 2, 2)
# eine einfache Methode, Zahlen zwischen 0 und 1 größer zu machen
# ist, die Quadratwurzel zu ziehen
ax.matshow(np.sqrt(the_picture_normalized)) 

plt.show()

Vektoren, Matrizen und Vektorisierung in Python 10

In [44]:

# Einen Blur-Filter über das Bild laufen lassen
# Das braucht etwas mehr Vorbereitung

# zunächst definieren wir den Filter, eine kleine Matrix
# dieser Filter mittelt n**2 benachbarte Werte in einem nxn-Pixel Quadrat
# das ist ein simpler Filter für Unschärfe
blur_size = 7
a_filter = 1/blur_size**2 * np.ones((blur_size, blur_size))

# nun müssen wir den Filter von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(the_picture_normalized)

# lege Array mit Nullen für das Resultat an
the_picture_normalized_blurred = np.zeros((y_dim, x_dim, 3))

# bestimme Abstand zum Rand
edge_dist = (blur_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # Loop über alle Kanäle
        for channel_index in np.arange(3):
            
            the_picture_normalized_blurred[y_index, x_index, channel_index] = np.sum( 
             a_filter * the_picture_normalized[y_index-edge_dist:y_index+edge_dist+1, 
                                                 x_index-edge_dist:x_index+edge_dist+1, 
                                                 channel_index])
    

# Das Bild mit Unschärfe darstellen
plt.matshow(the_picture_normalized_blurred)

plt.show()

100%|██████████| 1506/1506 [01:08<00:00, 21.84it/s]

Vektoren, Matrizen und Vektorisierung in Python 11

In [45]:

# Die Unschärfe ist im Ganzen schwer zu sehen, daher hier ein Ausschnitt

# plt.matshow(the_picture_normalized[300:700,700:1150,:])        # original
plt.matshow(the_picture_normalized_blurred[300:700,700:1150,:])  # blurred

plt.show()

In [46]:

# Einen allgemeinen Filter über das Bild laufen lassen
# Das braucht etwas mehr Vorbereitung

# zunächst definieren wir den Filter, eine kleine Matrix
# dieser Filter addiert benachbarte Werte mit bestimmten Gewichten 
# in einem nxn-Pixel Quadrat
# das ist ein simpler Filter für Kanten-Erkennung
filter_size = 3
a_filter = np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]])

# nun müssen wir den Filter von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(the_picture_normalized)

# lege Array mit Nullen für das Resultat an
the_picture_normalized_filtered = np.zeros((y_dim, x_dim, 3))

# bestimme Abstand zum Rand
edge_dist = (filter_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # Loop über alle Kanäle
        for channel_index in np.arange(3):
            
            the_picture_normalized_filtered[y_index, x_index, channel_index] = np.sum( 
             a_filter * the_picture_normalized[y_index-edge_dist:y_index+edge_dist+1, 
                                                 x_index-edge_dist:x_index+edge_dist+1, 
                                                 channel_index])
    

# Nochmal normieren, um sicher zu gehen, dass die Werte zwischen 0 und 1 bleiben
the_min = np.amin(the_picture_normalized_filtered)
the_max = np.amax(the_picture_normalized_filtered)
the_picture_normalized_filtered = (the_picture_normalized_filtered - the_min)/(the_max - the_min)


# Das gefilterte Bild darstellen
plt.matshow(the_picture_normalized_filtered)

plt.show()

100%|██████████| 1510/1510 [01:10<00:00, 21.30it/s]

Vektoren, Matrizen und Vektorisierung in Python 13

In [47]:

# Die Edge-Detection ist im Ganzen schwer zu sehen, daher hier ein Ausschnitt
# Außerdem eingeschränkt auf Kanal 2 (blau)

plt.matshow(the_picture_normalized_filtered[300:700,700:1150,2], cmap="Greys")

plt.show()

3.5 Was macht (und bringt) Vektorisierung?¶

Vektorisierung ist die Technik (und manchmal Kunst), Operationen, die man intuitiv als Loop ausführen würde, für Vektoren und Matrizen zu schreiben. Schafft man das, hat das im Allgemeinen Vorteile hinsichtlich der Geschwindigkeit. Der Vorteil endet allerdings, wenn der Speicher für die zu speichernden Arrays zu Ende geht.

Hier ein paar Beispiele.

In [48]:

# Die Anzahl der Punkte/die Dimension der Vektoren und Matrizen, mit der wir testen
n_points = 2000

# Erzeuge einen Vektor mit Zufallszahlen zwischen 0 und 1
a_random_vector = np.random.random(n_points)

In [49]:

%%time

# getimter Run: Quadriere alle Elemente im Vektor über einen Loop
a_random_vector_squared_1 = [elem**2 for elem in a_random_vector]

CPU times: user 710 µs, sys: 31 µs, total: 741 µs
Wall time: 751 µs

In [50]:

%%time

# getimter Run: Vektorisiertes elementweises Quadrat
a_random_vector_squared_2 = a_random_vector**2

CPU times: user 94 µs, sys: 51 µs, total: 145 µs
Wall time: 98 µs

In [51]:

# Erzeuge eine nxn Matrix mit Zufallszahlen zwischen 0 und 1
a_random_matrix = np.random.random((n_points, n_points))

In [52]:

%%time

# Initialisiere die Resultatmatrix
a_random_matrix_twice_1 = np.zeros((n_points, n_points))

# getimter Run: Verdopple alle Elemente der Matrix über zwei Loops
for ind_1 in range(n_points):
    for ind_2 in range(n_points):
        
        # Elementweises Setzen der Werte
        a_random_matrix_twice_1[ind_1, ind_2] = 2 * a_random_matrix[ind_1, ind_2]

CPU times: user 1.87 s, sys: 14.2 ms, total: 1.88 s
Wall time: 1.88 s

In [53]:

%%time

# getimter Run: Vektorisiertes Verdoppeln
a_random_matrix_twice_2 = 2 * a_random_matrix

CPU times: user 4.04 ms, sys: 6.68 ms, total: 10.7 ms
Wall time: 9.21 ms

In [54]:

%%time

# Hier noch etwas komplizierteres: Elementweise Abfrage

# Initialisiere die Resultatmatrix
a_random_matrix_gthalf_1 = np.zeros((n_points, n_points))

# getimter Run: Gehe über alle Elemente der Matrix über zwei Loops
for ind_1 in range(n_points):
    for ind_2 in range(n_points):
        
        # Zuerst den Wert an dieser Stelle der Matrix ausrechnen
        value_here = a_random_matrix[ind_1, ind_2]
        
        # if-Abfrage, ob der Wert größer als 0.5 ist
        if value_here > 0.5:
            
            # ja, ist größer, setze Resultat auf 1
            a_random_matrix_gthalf_1[ind_1, ind_2] = 1.
            
        else:
            # nein, ist kleiner, setze Resultat auf 0
            a_random_matrix_gthalf_1[ind_1, ind_2] = 0.

CPU times: user 1.39 s, sys: 11.8 ms, total: 1.4 s
Wall time: 1.4 s

In [55]:

%%time

# getimter Run: Vektorisierte if-Abfrage
# np.where(Bedingung, Resultat bei True, Resultat bei False)
a_random_matrix_gthalf_2 = np.where(a_random_matrix > 0.5 * np.ones((n_points, n_points)), 
                                    np.ones((n_points, n_points)), 
                                    np.zeros((n_points, n_points))
                                   )

CPU times: user 15 ms, sys: 22 ms, total: 37 ms
Wall time: 45.3 ms

In [56]:

# hier noch der Vergleich, ob die beiden Resultate auch gleich sind
# dazu summieren wir alle Stellen, wo es ungleich ist (also eine Eins/True steht)
# Die Summe ist also die Anzahl der falschen Vergleich-Elemente
np.sum(a_random_matrix_gthalf_1 != a_random_matrix_gthalf_2)

Out[56]:

Wichtig beim Umgang mit Arrays ist auch das Konzept der axis bei einer NumPy-Operation. Zum Beispiel können Sie die Summe aus einem Array komplett berechnen, aber auch nur entlang einer bestimmten Richtung (axis). Diese Vorgehensweise funktioniert bei den meisten Befehlen, die auf Arrays angewendet werden.

Hier auch dazu ein paar Beispiele:

In [57]:

# Nochmal die Matrix von oben:
a_matrix = np.array([[2, 3, 5, -3], [1, 3, -1, 4], [1, 2, -3, 2], [3, 4, -2, 1]])
a_matrix

Out[57]:

array([[ 2,  3,  5, -3],
       [ 1,  3, -1,  4],
       [ 1,  2, -3,  2],
       [ 3,  4, -2,  1]])

In [58]:

# Hier drei Summen: zunächst über alle Elemente
print(np.sum(a_matrix))

# Dann entlang der beiden Achsen separat
# Dabei entsteht ein Vektor mit den Summen entlang, z.B. der Zeilen
print(np.sum(a_matrix, axis=0))

# oder der Spalten
print(np.sum(a_matrix, axis=1))

# Dabei sollten die Summen auch zusammenpassen:
print(np.sum(np.sum(a_matrix, axis=1)))
print(np.sum(np.sum(a_matrix, axis=0)))

22
[ 7 12 -1  4]
[7 7 2 6]
22
22

In [59]:

# Erhöhen wir nun die Dimension ein bisschen
cubicle = np.concatenate((np.expand_dims(a_matrix, axis=0), np.expand_dims(np.transpose(a_matrix), axis=0)))
cubicle

Out[59]:

array([[[ 2,  3,  5, -3],
        [ 1,  3, -1,  4],
        [ 1,  2, -3,  2],
        [ 3,  4, -2,  1]],

       [[ 2,  1,  1,  3],
        [ 3,  3,  2,  4],
        [ 5, -1, -3, -2],
        [-3,  4,  2,  1]]])

In [60]:

# Hier kann man jetzt z.B. über zwei von drei Dimensionen summieren
# dazu muss man dem _axis_-Argument auch ein Tupel übergeben:
np.sum(cubicle, axis=(0, 2))

Out[60]:

array([14, 19,  1, 10])

3.6 Übungsaufgabe: Bildanalyse und einfache -manipulation¶

In dieser Aufgabe verwenden Sie bitte das Bild quadrate.png von der Moodle-Seite des Kurses. Die folgenden Aufgaben können Sie mit Hilfe der Befehle aus dieser Einheit gut erledigen:

Laden Sie das Bild in ein NumPy-Array. Hinweis: ein png hat 4 Kanäle: R, G, B, alpha. Den alpha-Kanal können Sie gleich nach dem Laden entfernen.
Zählen Sie die roten Quadrate (mit Python, nicht mit der Hand 🙂 ). Hinweis: Jedes Quadrat ist genau 81 Pixel groß.
Kopieren Sie das Array und färben Sie in der Kopie die roten Quadrate auf blau um

Zusatzaufgabe (optional): Bestimmen Sie die Mittelpunkte der Quadrate als Paare $i_y,i_x$ von Pixelpositionen des jeweiligen Quadrat-Mittelpunkts im Bild und geben Sie die Liste im Notebook aus.

In [61]:

# lade ein Bild über matplotlib-Bibliothek "image"
q_picture = image.imread(os.path.join('data', 'quadrate.png'))

# reduziere die 4 Kanäle aus dem png auf die üblichen 3
q_picture = q_picture[:,:,:-1]

# das Bild über matshow darstellen
plt.matshow(q_picture)

plt.show()

Vektoren, Matrizen und Vektorisierung in Python 15

In [62]:

# Die Lösung für die Aufgabe des Zählens der Quadrate kann man verschieden angehen
# die Basis dafür ist, die roten Pixel im Bild zu zählen und deren
# Anzahl dann durch 81 zu dividieren (die Größe eines Quadrats)

# um zwischen roten und weißen Pixeln zu unterscheiden, ist es wichtig, zu
# verstehen, dass weiße Pixel im bild die Kanalwerte (1,1,1) haben und
# einfärbige z.B. den Wert (1,0,0) 
# Wenn Sie also die Summe über die Pixel innerhalb eines Kanals nehmen,
# dann bekommen Sie die Summe aller Einsen. 
# in unserem Fall muss z.B. die Summer der Einsen im roten Kanal genau die 
# Gesamtzahl der Bildpixel ergeben.
# Die Summe der Einsen im grünen Kanal dagegen sollte um die Anzahl aller roten
# Bildpunkte geringer sein (weil dort im Kanal G eine Null statt einer Eins steht)

# Daher ist die Folgende Summe interessant: Über alle Dimensionen des Bildes, außer dem Kanal
channel_sums = np.sum(q_picture, axis=(0,1), dtype=int)
channel_sums

Out[62]:

array([2073600, 2071575, 2071575])

In [63]:

# damit erhalten wir nun die Lösung auf die Frage, wie viele rote Quadrate 
# es im Bild gibt, über den Folgenden Ausdruck: (Anzahl rot minus Anzahl grün)/81
(channel_sums[0] - channel_sums[1]) / 81

Out[63]:

25.0

In [64]:

# nun bleibt noch die zweite Aufgabe, nämlich die roten Punkte blau zu machen
# wir wissen bereits, dass bei roten Punkten die folgenden 
# Kanalwerte vorliegen: (1,0,0)
# Bei einem blauen Punkt hat man hingegen (0,0,1)
# Hier handelt es sich um eine zyklische Permutation, d.h. die Eins
# wird einfach zyklisch durchgetauscht.
# Da bei weißen Punkten (1,1,1) zyklische Vertauschungen keine Veränderung
# hervorrufen, kann man eine zyklische Vertauschung der kompletten
# Kanal-Ebenen im Bild verwenden, um die Farbe von rot auf blau zu ändern
# numpy hat für zyklische Vertauschungen in einem Array den Befehl np.roll bereitgestellt

# tausche die channels zyklisch, und zwar zu grün und blau, der Vollständigkeit halber
# Achten Sie auf axis=2, wodurch die Vertauschung auf die Kanaldimension beschränkt wird
green_q = np.roll(q_picture, 1, axis=2)
blue_q = np.roll(q_picture, 2, axis=2)

# jetzt plotten wir noch beide umgefärbten Bilder nebeneinander
fig = plt.figure(figsize=(10,5))

# die grünen Quadrate
ax = plt.subplot(1, 2, 1)
ax.matshow(green_q) 

# und die blauen Quadrate
ax = plt.subplot(1, 2, 2)
ax.matshow(blue_q) 


plt.show()

Vektoren, Matrizen und Vektorisierung in Python 16

In [65]:

# nun zur Zusatzaufgabe: wir möchten die Positionen der 
# Quadrate herauslesen. Dazu können wir im Prinzip das Filter-Prinzip und
# Teile des Codes von oben verwenden.
# Wir lassen einen Summenfilter der Größe 9x9 über z.B. den grünen Kanal laufen 
# und überall, wo er Null ergibt, entspricht die Filterposition dem Mittelpunkt des Quadrats

# den Filter brauchen wir hier nicht extra zu definieren, denn
# wir summieren einfach die Positionen im Farbkanal auf.
# Wir behalten aber das praktische Prinzip der Filtergröße bei
filter_size = 9


# wir müssen den Filter wieder von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(q_picture)

# bestimme Abstand zum Rand
edge_dist = (filter_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# lege Liste für Positionen an
position_list = []

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # aus dem Loop über alle Kanäle wird die Einschränkung auf den grünen (1)
        # wir wollen diesmal allerdings kein neues Bild erzeugen, sondern nur
        # die Summe checken
        this_sum = np.sum(q_picture[y_index-edge_dist:y_index+edge_dist+1, 
                                    x_index-edge_dist:x_index+edge_dist+1, 
                                    1])
        # checke, ob die Summe 0 ist
        if this_sum == 0:
            # ja, hänge die Position an die Liste an
            position_list.append([y_index, x_index])

# Hier die fertige Liste
print("Liste der Positionen: ", position_list)

# Checke die Länge der Liste (sollte 25 sein)
print("Länge dieser Liste: ", len(position_list))

100%|██████████| 1072/1072 [00:14<00:00, 73.15it/s]

Liste der Positionen:  [[5, 60], [16, 511], [27, 1534], [38, 1677], [60, 731], [170, 1006], [192, 698], [203, 1083], [280, 1028], [346, 456], [445, 1116], [467, 929], [478, 687], [500, 1699], [544, 445], [643, 1589], [654, 588], [698, 5], [720, 280], [797, 1578], [841, 401], [874, 324], [907, 918], [918, 764], [1072, 1523]]
Länge dieser Liste:  25

In [ ]:

Andreas Krassnigg: Physics, AI and Data Science

Theoretical Physicist at the University of Graz and Head of AI and Data Science at Innophore, Graz, Austria

Vektoren, Matrizen und Vektorisierung in Python