行列と連立一次方程式

ベクトルに作用する数の表

行列は数の長方形の表です。Python ではその行のリストとして持ち、[[1, 2], [3, 4]] は2×2行列、m[行][列] で1つの要素を読みます。行列は、同じ種類の計算を多数の行で一度に表します。各行がポートフォリオのウェイトで、ベクトルがその日の各資産のリターンなら、行列とベクトルの積は、すべてのポートフォリオのリターンを一度に与えます。結果の各要素は、ある行とベクトルの内積なので、行の長さとベクトルの長さは等しくなければなりません。

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)

連立一次方程式を解く

多くの問題は「A @ x = b となる x を求める」という形をしています。2つのヘッジ商品をそれぞれどれだけ持てば、あるエクスポージャーを打ち消せるか。どのファクターローディングなら観測されたリターンを再現できるか。これが連立一次方程式です。ガウスの消去法は、方程式を三角形に変えて解きます。各列で、その列の絶対値が最大の行をピボットに選び（上に入れ替えると丸め誤差が小さくなります）、その倍数を下の行から引いて0を作り、後退代入で下から上へ解きます。ピボットが（ほぼ）0なら行が従属していて行列は特異となり、答えは一意に決まりません。0で割らずにエラーを送出しましょう。

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 はピボット付き消去と後退代入で解く。ピボットが（ほぼ）0なら特異

次は確率変数と期待値です。あらゆるリスク・価格モデルの背後にある不確実性の言語です。

ラボ：採点される内容

• mat_vec で行列をベクトルに作用させる25点
• transpose で行と列を入れ替える25点
• mat_mul で2つの行列を掛ける25点
• solve で A x = b を消去法で解く25点