3. Lineare Algebra mit NumPy¶

3.1. NumPy installieren¶

In vielen mathematischen Anwendungen muss mit Vektoren und Matrizen gerechnet werden, beispielsweise bei der numerischen Berechnung von Integralen, Differentialgleichungen oder Problemen aus der Graphentheorie. Wir wollen uns in diesem Kapitel mit einer Python-Bibliothek beschäftigen, welche entsprechende Datentypen für Vektoren und Matrizen, sowie die üblichen Rechenoperationen für diese bereitstellt. Die Rede ist vom Paket NumPy.

Zunächst muss die NumPy-Bibliothek installiert werden. Je nachdem, welches System man verwendet, und abhängig vom persönlichen Geschmack verwendet man eine der 3 Varianten:

über den Paketmanager der Linux-Distribution (apt bzw. aptitude bei Ubuntu und Debian, zypper bei openSuse, yum bei Read Hat, pacman bei ArchLinux)

sudo apt-get install python3-numpy

über das Python-Paketverwaltungsprogramm pip

pip3 install -U numpy

über die Conda-Umgebung

conda install numpy

Nun muss die Bibliothek in unser Python-Skript eingebunden werden. Wir könnten dies wie in Abschnitt Variablen und Datentypen mit der Zeile from numpy import * machen, was aber nicht empfehlenswert ist, da die NumPy-Bibliothek die Funktionen sin, cos, sqrt, etc. bereit stellt, welche die aus der math-Bibliothek überschreiben würden. Daher nutzen wir folgende Variante:

import numpy as np

Der Zusatz as np gibt nur an, dass wir die Bibliothek in Zukunft unter dem kürzeren Namen np und nicht unter dem langen numpy ansprechen können.

3.2. Arbeiten mit Vektoren¶

3.2.1. Vektoren erzeugen¶

Ein Vektor bzw. NumPy-Array kann mit der Funktion np.array(...) von einem beliebigen iterierbaren Objekt (Liste, Tupel, …), deren Elemente vom gleichen Typ sind, initialisiert werden:

a = np.array([1.,2.,3.])
b = np.array((6.,5.,4.))
print("a ist ein", type(a), "mit Wert", a)
print("b ist ein", type(b), "mit Wert", b)

a ist ein <class 'numpy.ndarray'> mit Wert [1. 2. 3.]
b ist ein <class 'numpy.ndarray'> mit Wert [6. 5. 4.]

Die Klasse ndarray repräsentiert ein mehrdimensionales Array. In unserem Fall ein Vektor, also ein Array der Dimension

a.ndim

mit Daten vom Typ

a.dtype

dtype('float64')

der Dimension

a.shape

(3,)

2 weitere Funktionen um NumPy-Arrays zu erzeugen sind:

x = np.linspace(0,3,11) # Äquidistantes Punktgitter für [0,3] aus 10 Punkten
x

array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3. ])

x = np.arange(0, 3, 0.3) # Äquidistantes Punktgitter für [0,3) mit Inkrement 0.25
x

array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7])

3.2.2. Elementare Vektoroperationen¶

Ein ndarray ist, wie eine Liste oder ein Tupel, ein iterierbares Objekt. Wir können also einfach mit einer for-Schleife über alle Elemente gehen:

val = 0
for e in a:
    val += e
print("Die Summe der Einträge von", a, "ist", val)

Die Summe der Einträge von [1. 2. 3.] ist 6.0

Wir können einige elementare Rechenoperationen für mit unseren Vektoren durchführen und es kommt das erwartete Ergebnis raus:

a+b

array([7., 7., 7.])

b-a

array([5., 3., 1.])

a*b

array([ 6., 10., 12.])

a/b

array([0.16666667, 0.4       , 0.75      ])

Wir beobachten, dass die Grundrechenarten einfach elementweise durchgeführt werden. Bei Addition und Subtraktion ist es auch das, was wir aus der Linearen Algebra-Perspektive erwarten würden, aber bei der Multiplikation und Division ist dies nicht klar. Eventuell hätten wir die Berechnung des Skalar- oder Kreuzproduktes erwartet. Dazu später mehr.

Die Rechenoperation +, -, *, / werfen auch einen Fehler, wenn die Größen der beiden Vektoren nicht kompatibel sind:

c = np.array([8,7,6,5])
a+c

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [13], in <module>
      1 c = np.array([8,7,6,5])
----> 2 a+c

ValueError: operands could not be broadcast together with shapes (3,) (4,) 

Tippen wir a.<TAB> ein um zu erfahren, welche Funktionen die Klasse ndarray bereitstellt. Hier einige Beispiele:

print("Der kleinste Eintrag von a ist", a.min(), "bei Index", a.argmin())
print("Das Skalarprodukt <a,b> ist", a.dot(b))
print("Die Summe aller Einträge von a ist", a.sum())

Der kleinste Eintrag von a ist 1.0 bei Index 0
Das Skalarprodukt <a,b> ist 28.0
Die Summe aller Einträge von a ist 6.0

3.2.3. Weitere Rechenoperationen¶

Neben diesen elementaren Funktionen, welche die Klasse ndarray direkt bereitstellt, finden wir noch weitere Rechenoperation, welche als freie Funktionen in der Bibliothek numpy implementiert sind. Wir tippen np.<TAB> ein und erhalten eine Liste all dieser Funktionen. Uns fällt auf, dass hier nochmals die Funktionen sqrt, exp, sin, cos definiert sind. Diese sind zwar schon in der math-Bibliothek vorhanden, lassen sich aber nicht auf Objekte vom Typ ndarray anwenden:

import math
math.exp(a)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [15], in <module>
      1 import math
----> 2 math.exp(a)

TypeError: only size-1 arrays can be converted to Python scalars

Die Exponentialfunktion aus der numpy-Bibliothek hingegen wendet die Exponentialfunktion komponentenweise auf den Vektor an:

np.exp(a)

array([ 2.71828183,  7.3890561 , 20.08553692])

Auch die Funktionen log, sin, cos, tan, abs, …, werden elementweise auf den Vektor angewendet. Genau so werden auch Vergleichsoperationen elementweise angewendet:

x = np.array([0,1,2,3,4])
y = (x <=2)
y

array([ True,  True,  True, False, False])

Übungsaufgabe

Berechne die Werte von \(\sin(x)\) für \(x\in \{0,\frac\pi6,\frac\pi3,\frac\pi2,\ldots,2\pi\}\).

Neben diesen elementaren Rechenoperationen finden wir im Modul numpy aber auch vektorspezifische Operationen, wie beispielsweise das Skalar- und Kreuzprodukt:

np.dot(a,b)

28.0

np.cross(a,b)

array([-7., 14., -7.])

Bei einer genaueren Betrachtung der Funktionen aus np stellt man fest, dass Funktionen für beispielsweise Vektornormen fehlen. Diese finden wir im Submodul numpy.linalg. Wir tippen np.linalg.<TAB> ein und erhalten wieder eine Liste aller angebotenen Funktionen.

Wir können beispielsweise wie folgt die üblichen Normen berechnen:

print("Euklidische Norm :", np.linalg.norm(a))
print("Maximumnorm      :", np.linalg.norm(a, np.Inf))
print("1-Norm           :", np.linalg.norm(a, 1))

Euklidische Norm : 3.7416573867739413
Maximumnorm      : 3.0
1-Norm           : 6.0

Vektoroperationen

Zusammenfassend stellen wir fest, dass die NumPy-Bibliothek sehr viele Rechenoperationen für Vektoren aus der linearen Algebra bereitstellt. Diese sind entweder

Klassenfunktionen von ndarray (z.B. a.min())
Freie Funktionen im Paket numpy (z.B. np.dot(a,b))
Freie Funktionen im Paket numpy.linalg (z.B. np.linalg.norm(a))

Übungsaufgabe

Schreibe eine Funktion, die den Winkel zweier Vektoren über die Formel

\[ \cos(\alpha) = \frac{u^\top v}{\|u\|\,\|v\|} \]

berechnet.

3.2.4. Zugriff auf Vektoreinträge¶

Schauen wir uns zuletzt noch an, wie auf einzelne oder mehrere Elemente des Vektors in einem bestimmten Bereich zugegriffen werden kann. Dies geschieht mit dem []-Operator:

a = np.linspace(0,1,11)
a

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

print("Vierter Eintrag:") # Die Zählung des Index' beginnt bei 0
a[3]

Vierter Eintrag:

0.30000000000000004

Wir können auch einen Teil des Arrays extrahieren, indem wir einen Indexbereich i:j angeben:

a[3:6] # Indizes zwischen 3 und 5

array([0.3, 0.4, 0.5])

a[:6] # Indizes bis 5

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5])

a[6:] # Indizes ab 6

array([0.6, 0.7, 0.8, 0.9, 1. ])

Wir können auch Teilvektoren überschreiben, natürlich unter Beachtung der Dimension:

a[1:4] = np.ones((3,))  # Schreibe Einsen in die Einträge 1 bis 3
a[6:8] = np.zeros((2,)) # Schreibe Nullen in die Einträge 6 und 7
a

array([0. , 1. , 1. , 1. , 0.4, 0.5, 0. , 0. , 0.8, 0.9, 1. ])

Nebenbei haben wir hier auch 2 Methoden kennengelernt um Arrays mit Einträgen 0 oder 1 zu initialisieren.

3.3. Arbeiten mit Matrizen¶

3.3.1. Erzeugen von Matrizen¶

Auch für Matrizen nutzen wir den Datentyp numpy.ndarray, mit dem Unterschied, dass wir ein zweidimensionales Array verwenden. Wir können Matrizen aus zweidimensionalen Listen (also Listen deren Einträge selbst wieder Listen sind) erzeugen:

A = np.array([[1.,4.,2.],[2.,0.,3.],[1.,1.,4.]])
A

array([[1., 4., 2.],
       [2., 0., 3.],
       [1., 1., 4.]])

Um eine Einheitsmatrix zu erzeugen nutzen wir

B = np.eye(3)
B

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Wir können auch eine leere Matrix erzeugen und die Einträge anschließend direkt eintragen:

C = np.zeros((3,3)) # Das Argument ist ein Tupel und gibt die Größer der Matrix an
C[0,0] = 1
C[1,1] = 2
C[1,2] = 5
C[2,1] = 4
C[2,2] = 1
C

array([[1., 0., 0.],
       [0., 2., 5.],
       [0., 4., 1.]])

Zu einem gegebenen Vektor kann man auch die zugehörige Diagonalmatrix erzeugen mit

D = np.diag(np.array([1.,2.,3.]))
D

array([[1., 0., 0.],
       [0., 2., 0.],
       [0., 0., 3.]])

3.3.2. Rechnen mit Matrizen¶

Wie schon bei den Vektoren beobachtet werden die Grundrechenarten sowie Funktionen exp, sin, cos, log, etx. aus dem Paket numpy elementweise angewendet:

C = A+B
C

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

D = A*B
D

array([[1., 0., 0.],
       [0., 0., 0.],
       [0., 0., 4.]])

E = np.exp(A)
E

array([[ 2.71828183, 54.59815003,  7.3890561 ],
       [ 7.3890561 ,  1.        , 20.08553692],
       [ 2.71828183,  2.71828183, 54.59815003]])

Weitere Rechenoperationen werden von der Klasse numpy.ndarray bereitgestellt, wie das Transponieren einer Matrix:

C.transpose()

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

Die Matrixmultiplikation wird als freie Funktion vom Modul numpy bereitgestellt:

np.matmul(A,B)

array([[1., 4., 2.],
       [2., 0., 3.],
       [1., 1., 4.]])

Auch in der Bibliothek np.linalg finden wir weitere Funktionen aus der linearen Algebra, beispielsweise eine Funktion zur Berechnung der Determinante

np.linalg.det(C)

-22.000000000000004

und der Inversen

np.linalg.inv(C)

array([[-0.09090909,  0.81818182, -0.45454545],
       [ 0.31818182, -0.36363636,  0.09090909],
       [-0.04545455, -0.09090909,  0.27272727]])

Die Eigenwerte- und Vektoren berechnen wir mit

E,V = np.linalg.eig(C)
print("Eigenwerte    :\n", E)
print("Eigenvektoren :\n", V)

Eigenwerte    :
 [ 6.97413644 -1.33574651  2.36161007]
Eigenvektoren :
 [[ 0.63932126  0.77319109 -0.8514755 ]
 [ 0.50515119 -0.63379132 -0.29406667]
 [ 0.57973321 -0.02200211  0.43418229]]

Die Eigenvektoren stehen spaltenweise in V. Testen wir die Rechnung indem wir \(C\,v - \lambda\,v=0\) überprüfen:

for i in range(3):
    res = np.matmul(C, V[:,i]) - E[i]*V[:,i]  # C*v-lambda*v    
    print("Fehler: ", np.linalg.norm(res))

Fehler:  3.66205343881779e-15
Fehler:  1.5852031831917945e-15
Fehler:  2.434909185329438e-15

Lineare Gleichungssysteme

numpy.linalg stellt auch Methoden zur Lösung linearer Gleichungssysteme bereit:

b = np.array([24.,23.,30.])
x = np.linalg.solve(C, b)
x

array([3., 2., 5.])

Eine einfache Probe ergibt

np.matmul(C,x)-b

array([ 0.00000000e+00,  0.00000000e+00, -3.55271368e-15])

Zugriff auf Matrixeinträge

Auch der Zugriff auf einzelne Einträge bzw. auf Teilmatrizen funktioniert analog zu eindimensionalen NumPy-Arrays:

print("Eintrag C_11           :", C[0,0])
print("Zweite Spalte          :", C[:,1])
print("Dritte Zeile           :", C[2,:])
print("Letzte (=dritte) Zeile :", C[-1,:])

Eintrag C_11           : 2.0
Zweite Spalte          : [4. 1. 1.]
Dritte Zeile           : [1. 1. 5.]
Letzte (=dritte) Zeile : [1. 1. 5.]

Übungsaufgabe

Schreibe eine Funktion, welche überprüft, ob eine gegebene Matrix eine Nullzeile besitzt.

Matrizen stapeln

Weiterhin lassen sich mit den Befehlen numpy.vstack und numpy.hstack Matrizen “zusammenkleben”:

np.vstack([A,B,C]) # Stapelt A, B, C vertikal

array([[1., 4., 2.],
       [2., 0., 3.],
       [1., 1., 4.],
       [1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.],
       [2., 4., 2.],
       [2., 1., 3.],
       [1., 1., 5.]])

np.hstack([A,B,C]) # Stapelt A, B, C vertikal

array([[1., 4., 2., 1., 0., 0., 2., 4., 2.],
       [2., 0., 3., 0., 1., 0., 2., 1., 3.],
       [1., 1., 4., 0., 0., 1., 1., 1., 5.]])

Um einen Vektor an eine Matrix zu kleben müssen wir bei numpy.hstack den Vektor zunächst in eine \(n\times 1\)-Matrix umwandeln wir folgendes Beispiel zeigt:

np.vstack([A,b])

array([[ 1.,  4.,  2.],
       [ 2.,  0.,  3.],
       [ 1.,  1.,  4.],
       [24., 23., 30.]])

np.hstack([A,b.reshape((3,1))])

array([[ 1.,  4.,  2., 24.],
       [ 2.,  0.,  3., 23.],
       [ 1.,  1.,  4., 30.]])

Mutable oder Immutable?

Zuletzt wollen wir noch überprüfen, wie sich NumPy-Arrays als Funktionsparameter verhalten. Sind diese mutable oder immutable? Ein einfacher Test gibt:

def modify_matrix(A):
    A[1,1] = 1.
    
A = np.zeros((3,3))
modify_matrix(A)
A

array([[0., 0., 0.],
       [0., 1., 0.],
       [0., 0., 0.]])

Objekte vom Typ numpy.ndarray sind offensichtlich mutable. Man muss an dieser Stelle auch bei einer Zuweisung von Matrizen aufpassen:

B = A
B[1,1] = 2.
A

array([[0., 0., 0.],
       [0., 2., 0.],
       [0., 0., 0.]])

Wir haben hier den Namen B an das gleiche Objekt wie A gebunden. Nachdem wir B verändert haben, hat sich diese Änderung offensichtlich auch auf die Matrix A ausgewirkt. Es wird also keine Kopie der Matrix angezeigt.

Möchte man tatsächlich eine Kopie einer Matrix erstellen, so nutzt man

B = np.copy(A)
B[1,1] = 3.
A

array([[0., 0., 0.],
       [0., 2., 0.],
       [0., 0., 0.]])

Wir haben hier B verändert, die Änderung wirkt sich aber nicht auf die Matrix A aus.

Übungsaufgabe

Erstelle die Matrix

\[\begin{split} A = \begin{pmatrix} 1 & & & & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & & 1 \end{pmatrix} \end{split}\]

welche bei der Finite-Differenzen-Diskretisierung des Randwertproblems \(-y''(t) = f(t)\) für \(t\in(0,1)\) und \(y(0)=y(1)=0\) auftritt.

Mathematisches Programmieren

Lineare Algebra mit NumPy

Contents

3. Lineare Algebra mit NumPy¶

3.1. NumPy installieren¶

3.2. Arbeiten mit Vektoren¶

3.2.1. Vektoren erzeugen¶

3.2.2. Elementare Vektoroperationen¶

3.2.3. Weitere Rechenoperationen¶

3.2.4. Zugriff auf Vektoreinträge¶

3.3. Arbeiten mit Matrizen¶

3.3.1. Erzeugen von Matrizen¶

3.3.2. Rechnen mit Matrizen¶