Matrizen und lineare Gleichungssysteme

Eine Zahlentabelle, die auf Vektoren wirkt

Eine Matrix ist eine rechteckige Zahlentabelle. In Python speichern wir sie als Liste von Zeilen, [[1, 2], [3, 4]] ist also eine 2×2-Matrix und m[zeile][spalte] liest einen Eintrag. Matrizen beschreiben viele Zeilen derselben Art von Rechnung auf einmal. Enthält jede Zeile die Gewichte eines Portfolios und der Vektor die Tagesrenditen der Werte, ergibt das Produkt aus Matrix und Vektor die Rendite jedes Portfolios in einem Schritt. Jeder Eintrag des Ergebnisses ist das Skalarprodukt einer Zeile mit dem Vektor, die Zeilenlänge muss also der Vektorlänge entsprechen.

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def mat_vec(matrix, vector):
    return [dot(row, vector) for row in matrix]

print(mat_vec([[1, 2], [3, 4]], [5, 6]))   # [17, 39]

# Two portfolios, three assets, one day of returns
weights = [[0.5, 0.3, 0.2],
           [0.2, 0.5, 0.3]]
returns = [0.02, -0.01, 0.03]
print([round(r, 4) for r in mat_vec(weights, returns)])   # [0.013, 0.008]

Transponierte und Matrixmultiplikation

Die Transponierte vertauscht Zeilen und Spalten: Eintrag (i, j) wird zu (j, i), aus einer 2×3-Matrix wird also 3×2. In Python erledigt das [list(spalte) for spalte in zip(*m)]. Das Produkt A @ B einer n×k- und einer k×m-Matrix ist eine n×m-Matrix, deren Eintrag (i, j) das Skalarprodukt von Zeile i von A mit Spalte j von B ist. Die Spalten von B sind die Zeilen ihrer Transponierten, was den Code kurz macht. Die inneren Größen müssen passen (k und k), und anders als bei Zahlen kommt es auf die Reihenfolge an: Im Allgemeinen ist A @ B nicht B @ A.

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def transpose(m):
    return [list(column) for column in zip(*m)]

def mat_mul(a, b):
    columns = transpose(b)
    return [[dot(row, column) for column in columns] for row in a]

print(transpose([[1, 2, 3], [4, 5, 6]]))   # [[1, 4], [2, 5], [3, 6]]

A = [[1, 2], [3, 4]]
B = [[5, 6], [7, 8]]
print(mat_mul(A, B))   # [[19, 22], [43, 50]]
print(mat_mul(B, A))   # [[23, 34], [31, 46]]  (a different answer: order matters)

Ein lineares Gleichungssystem lösen

Viele Fragen haben die Form: Finde x mit A @ x = b. Wie viel von jedem von zwei Absicherungsinstrumenten hebt ein gegebenes Exposure auf? Welche Faktorladungen reproduzieren beobachtete Renditen? Das ist ein lineares Gleichungssystem. Das Gauß-Verfahren löst es, indem es die Gleichungen in ein Dreieck bringt: Für jede Spalte wählt man die Zeile mit dem größten Betrag in dieser Spalte als Pivot (nach oben getauscht, das hält Rundungsfehler klein), zieht Vielfache davon von den Zeilen darunter ab, um Nullen zu erzeugen, und löst dann von unten nach oben per Rückwärtseinsetzen. Ist ein Pivot (fast) null, sind die Zeilen abhängig, die Matrix ist singulär und es gibt keine eindeutige Antwort: Lösen Sie einen Fehler aus, statt durch null zu teilen.

def solve(a, b):
    n = len(a)
    m = [list(map(float, row)) + [float(b[i])] for i, row in enumerate(a)]  # augmented matrix
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(m[r][col]))
        if abs(m[pivot][col]) < 1e-12:
            raise ValueError("the matrix is singular")
        m[col], m[pivot] = m[pivot], m[col]
        for r in range(col + 1, n):
            factor = m[r][col] / m[col][col]
            for c in range(col, n + 1):
                m[r][c] -= factor * m[col][c]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):          # back-substitution
        known = sum(m[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (m[i][n] - known) / m[i][i]
    return x

print(solve([[2, 1], [1, 3]], [3, 5]))                    # [0.8, 1.4]
print(solve([[1, 1, 1], [0, 2, 5], [2, 5, -1]], [6, -4, 27]))  # [5.0, 3.0, -2.0]

try:
    solve([[1, 2], [2, 4]], [1, 2])         # the second row is twice the first
except ValueError as error:
    print(error)                            # the matrix is singular

import math

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def mat_vec(matrix, vector):
    return [dot(row, vector) for row in matrix]

# Check a solution by multiplying back: A x should give b again
A = [[2, 1], [1, 3]]
b = [3, 5]
x = [0.8, 1.4]                       # the answer from solve(A, b)
residual = [got - want for got, want in zip(mat_vec(A, x), b)]
print([abs(round(r, 12)) for r in residual])  # [0.0, 0.0]  (tiny rounding errors vanish)
print(all(math.isclose(g, w) for g, w in zip(mat_vec(A, x), b)))   # True

# The identity matrix leaves any vector unchanged
identity = [[1, 0], [0, 1]]
print(mat_vec(identity, [7, 9]))     # [7, 9]

• Eine Matrix ist eine Liste von Zeilen; mat_vec bildet das Skalarprodukt jeder Zeile mit dem Vektor.
• Transponieren mit zip(*m); A @ B Eintrag (i, j) = Zeile i von A mal Spalte j von B; die Reihenfolge zählt.
• Lösen Sie A x = b per Elimination mit Pivot und Rückwärtseinsetzen; ein (fast) Null-Pivot heißt singulär.

Als Nächstes: Zufallsvariablen und Erwartungswert, die Sprache der Unsicherheit hinter jedem Risiko- und Bewertungsmodell.

Lab: Was bewertet wird

• Mit mat_vec eine Matrix auf einen Vektor anwenden25 Pkt.
• Mit transpose Zeilen und Spalten vertauschen25 Pkt.
• Mit mat_mul zwei Matrizen multiplizieren25 Pkt.
• Mit solve A x = b durch Elimination lösen25 Pkt.