A table of numbers that acts on vectors
A matrix is a rectangular table of numbers. In Python we store it as a list of rows, so [[1, 2], [3, 4]] is a 2×2 matrix and m[row][column] reads one entry. Matrices describe many rows of the same kind of calculation at once. If each row holds the weights of a portfolio and the vector holds the day's asset returns, the product of the matrix and the vector is the return of every portfolio in one step. Each result entry is the dot product of one row with the vector, so the row length must equal the vector length.
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]Transpose and matrix multiplication
The transpose swaps rows and columns: entry (i, j) becomes (j, i), so a 2×3 matrix becomes 3×2. In Python [list(column) for column in zip(*m)] does it. The product A @ B of an n×k and a k×m matrix is an n×m matrix whose entry (i, j) is the dot product of row i of A with column j of B. Columns of B are the rows of its transpose, which makes the code short. The inner sizes must match (k and k), and, unlike numbers, the order matters: in general A @ B is not 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)Solving a linear system
Many questions have the form: find x such that A @ x = b. How much of each of two hedge instruments cancels a given exposure? Which factor loadings reproduce observed returns? This is a linear system. Gaussian elimination solves it by turning the equations into a triangle: for each column, pick the row with the largest absolute value in that column as the pivot (swap it up, which keeps rounding errors small), subtract multiples of it from the rows below to make zeros, then solve upwards by back-substitution. If a pivot is (almost) zero the rows are dependent, the matrix is singular, and there is no single answer: raise an error instead of dividing by zero.
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 singularimport 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]x = A⁻¹ b, but in code you almost never compute the inverse: eliminating directly is faster and more accurate. Also check the answer by multiplying back (A @ x should reproduce b to rounding), and treat a nearly singular matrix, such as two assets that are almost perfectly correlated, as a warning that the solution will be unstable. Later you will use numpy.linalg.solve, which does the same work in compiled code.- A matrix is a list of rows;
mat_vectakes the dot product of every row with the vector. - Transpose with
zip(*m);A @ Bentry (i, j) = row i of A dot column j of B; the order matters. - Solve
A x = bby elimination with pivoting and back-substitution; a (near) zero pivot means singular.
Next up: random variables and expectation, the language of uncertainty behind every risk and pricing model.
Lab · what's graded
- Apply a matrix to a vector with
mat_vec25 pts - Swap rows and columns with
transpose25 pts - Multiply two matrices with
mat_mul25 pts - Solve A x = b by elimination with
solve25 pts