Matrices y sistemas lineales

Una tabla de números que actúa sobre vectores

Una matriz es una tabla rectangular de números. En Python la guardamos como una lista de filas, así que [[1, 2], [3, 4]] es una matriz 2×2 y m[fila][columna] lee una entrada. Las matrices describen muchas filas del mismo tipo de cálculo a la vez. Si cada fila contiene los pesos de una cartera y el vector los rendimientos del día de los activos, el producto de la matriz por el vector da el rendimiento de cada cartera en un paso. Cada entrada del resultado es el producto escalar de una fila con el vector, así que la longitud de la fila debe igualar la del vector.

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]

Transpuesta y multiplicación de matrices

La transpuesta intercambia filas y columnas: la entrada (i, j) pasa a (j, i), así que una matriz 2×3 pasa a 3×2. En Python lo hace [list(columna) for columna in zip(*m)]. El producto A @ B de una matriz n×k y una k×m es una matriz n×m cuya entrada (i, j) es el producto escalar de la fila i de A con la columna j de B. Las columnas de B son las filas de su transpuesta, lo que acorta el código. Las dimensiones interiores deben coincidir (k y k) y, a diferencia de los números, el orden importa: en general A @ B no es 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)

Resolver un sistema lineal

Muchas preguntas tienen la forma: hallar x tal que A @ x = b. ¿Cuánto de cada uno de dos instrumentos de cobertura anula una exposición dada? ¿Qué cargas de factores reproducen los rendimientos observados? Es un sistema lineal. La eliminación gaussiana lo resuelve convirtiendo las ecuaciones en un triángulo: para cada columna se elige como pivote la fila con mayor valor absoluto en esa columna (se sube, lo que mantiene pequeños los errores de redondeo), se restan múltiplos de ella a las filas inferiores para crear ceros y luego se resuelve hacia arriba por sustitución regresiva. Si un pivote es (casi) cero, las filas son dependientes, la matriz es singular y no hay una respuesta única: lanza un error en lugar de dividir por cero.

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]

• Una matriz es una lista de filas; mat_vec hace el producto escalar de cada fila con el vector.
• Transpuesta con zip(*m); A @ B entrada (i, j) = fila i de A por columna j de B; el orden importa.
• Resuelve A x = b por eliminación con pivote y sustitución regresiva; un pivote (casi) cero significa singular.

A continuación: variables aleatorias y esperanza, el lenguaje de la incertidumbre detrás de todo modelo de riesgo y valoración.

Lab: qué se califica

• Aplicar una matriz a un vector con mat_vec25 pts
• Intercambiar filas y columnas con transpose25 pts
• Multiplicar dos matrices con mat_mul25 pts
• Resolver A x = b por eliminación con solve25 pts