矩阵与线性方程组

作用于向量的数表

矩阵是由数字构成的矩形表格。在 Python 中我们把它存成行的列表，所以 [[1, 2], [3, 4]] 是 2×2 矩阵，m[行][列] 读取一个元素。矩阵可以一次性描述许多行同类的计算。如果每一行是一个投资组合的权重，向量是各资产当天的收益率，那么矩阵与向量的乘积就一步给出每个组合的收益。结果的每个元素都是某一行与向量的点积，所以行长必须等于向量长度。

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]

转置与矩阵乘法

转置交换行与列：元素 (i, j) 变成 (j, i)，因此 2×3 矩阵变成 3×2。在 Python 中，[list(column) for column in zip(*m)] 就能做到。n×k 矩阵与 k×m 矩阵的乘积 A @ B 是一个 n×m 矩阵，其元素 (i, j) 是 A 的第 i 行与 B 的第 j 列的点积。B 的列就是其转置的行，这让代码很短。内部维度必须相同（k 与 k），而且与数字不同，顺序很重要：一般来说 A @ B 不等于 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)

求解线性方程组

很多问题都是这样的形式：求 x 使 A @ x = b。两种对冲工具各要多少才能抵消给定的敞口？哪些因子载荷能复现观察到的收益？这就是线性方程组。高斯消元法把方程化成三角形来求解：对每一列，选出该列绝对值最大的行作为主元（换到上面，可使舍入误差较小），用它的倍数减去下面各行以产生零，再通过回代自下而上求解。如果主元（几乎）为零，说明各行线性相关，矩阵是奇异的，没有唯一解：应报错，而不是除以零。

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]

• 矩阵是行的列表；mat_vec 取每一行与向量的点积。
• 用 zip(*m) 转置；A @ B 的 (i, j) 元素 = A 的第 i 行点乘 B 的第 j 列；顺序有关。
• 用带主元的消元和回代解 A x = b；主元（接近）为零表示奇异。

下一节：随机变量与期望，所有风险与定价模型背后的不确定性语言。

实验：评分内容

• 用 mat_vec 把矩阵作用于向量25 分
• 用 transpose 交换行与列25 分
• 用 mat_mul 把两个矩阵相乘25 分
• 用 solve 以消元法求解 A x = b25 分