Matrices et systèmes linéaires

Un tableau de nombres qui agit sur des vecteurs

Une matrice est un tableau rectangulaire de nombres. En Python, on la stocke comme une liste de lignes : [[1, 2], [3, 4]] est une matrice 2×2 et m[ligne][colonne] lit une entrée. Les matrices décrivent d'un coup de nombreuses lignes d'un même calcul. Si chaque ligne contient les poids d'un portefeuille et que le vecteur contient les rendements du jour, le produit de la matrice par le vecteur donne le rendement de chaque portefeuille en une seule étape. Chaque entrée du résultat est le produit scalaire d'une ligne avec le vecteur, donc la longueur d'une ligne doit égaler celle du vecteur.

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]

Transposée et produit matriciel

La transposée échange lignes et colonnes : l'entrée (i, j) devient (j, i), donc une matrice 2×3 devient 3×2. En Python, [list(colonne) for colonne in zip(*m)] le fait. Le produit A @ B d'une matrice n×k et d'une matrice k×m est une matrice n×m dont l'entrée (i, j) est le produit scalaire de la ligne i de A avec la colonne j de B. Les colonnes de B sont les lignes de sa transposée, ce qui raccourcit le code. Les dimensions intérieures doivent coïncider (k et k) et, contrairement aux nombres, l'ordre compte : en général A @ B n'est pas 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)

Résoudre un système linéaire

Beaucoup de questions prennent la forme : trouver x tel que A @ x = b. Quelle quantité de chacun de deux instruments de couverture annule une exposition donnée ? Quelles expositions aux facteurs reproduisent les rendements observés ? C'est un système linéaire. L'élimination de Gauss le résout en transformant les équations en un triangle : pour chaque colonne, on choisit comme pivot la ligne ayant la plus grande valeur absolue dans cette colonne (on la remonte, ce qui limite les erreurs d'arrondi), on retranche des multiples de cette ligne aux lignes du dessous pour créer des zéros, puis on résout en remontant par substitution arrière. Si un pivot est (presque) nul, les lignes sont dépendantes, la matrice est singulière et il n'y a pas de réponse unique : levez une erreur plutôt que de diviser par zéro.

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]

• Une matrice est une liste de lignes ; mat_vec fait le produit scalaire de chaque ligne avec le vecteur.
• Transposée avec zip(*m) ; A @ B entrée (i, j) = ligne i de A scalaire colonne j de B ; l'ordre compte.
• Résolvez A x = b par élimination avec pivot et substitution arrière ; un pivot (presque) nul signifie singulière.

À suivre : variables aléatoires et espérance, le langage de l'incertitude derrière tout modèle de risque et de valorisation.

Lab : ce qui est noté

• Appliquer une matrice à un vecteur avec mat_vec25 pts
• Échanger lignes et colonnes avec transpose25 pts
• Multiplier deux matrices avec mat_mul25 pts
• Résoudre A x = b par élimination avec solve25 pts