Matrix Operations

Matrix Multiplication: The Rules

Matrix multiplication is not something you do entry-by-entry. It is a deliberate combination of rows and columns, and it has specific rules that differ from everyday arithmetic.

To multiply $\mathbf{A}$ (size $m \times n$) by $\mathbf{B}$ (size $n \times p$), you produce an $m \times p$ matrix. The entry at row $i$, column $j$ is the dot product of row $i$ from $\mathbf{A}$ with column $j$ from $\mathbf{B}$:

This has two immediate consequences. First, the inner dimensions must match, the number of columns in $\mathbf{A}$ must equal the number of rows in $\mathbf{B}$. Second, matrix multiplication is not commutative: $\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}$ in general. The order matters.

Matrix multiplication is associative: $(\mathbf{A}\mathbf{B})\mathbf{C} = \mathbf{A}(\mathbf{B}\mathbf{C})$, and distributive: $\mathbf{A}(\mathbf{B} + \mathbf{C}) = \mathbf{A}\mathbf{B} + \mathbf{A}\mathbf{C}$.

The Identity Matrix and Inverses

The identity matrix $\mathbf{I}$ plays the same role in matrix multiplication that 1 plays in scalar multiplication. It is a square matrix with 1s on the diagonal and 0s everywhere else:

For any compatible matrix $\mathbf{A}$: $\mathbf{A}\mathbf{I} = \mathbf{A}$ and $\mathbf{I}\mathbf{A} = \mathbf{A}$.

The inverse of a square matrix $\mathbf{A}$, written $\mathbf{A}^{-1}$, is the matrix that satisfies:

Not all matrices have inverses. A matrix that cannot be inverted is called singular or non-invertible. Singularity is not a rare edge case, it is a warning sign that your matrix contains redundant information, like two identical columns.

The determinant $\det(\mathbf{A})$ is a scalar derived from the matrix that tells you whether inversion is possible. If $\det(\mathbf{A}) = 0$, the matrix is singular. If $\det(\mathbf{A}) \neq 0$, an inverse exists.

Properties of Inverses and Transposes

A few algebraic rules are used constantly in derivations:

Notice the order reverses. To undo the operation $\mathbf{A}$ then $\mathbf{B}$, you first undo $\mathbf{B}$ then $\mathbf{A}$.

Again, order reverses when you transpose a product.

The trace of a matrix is the sum of its diagonal elements. It has the useful property of being invariant under transposition: $\text{tr}(\mathbf{A}) = \text{tr}(\mathbf{A}^\top)$.

Solving Linear Systems

One of the most common problems in quantitative finance is solving for unknowns in a system of linear equations. This problem can be written compactly as:

where $\mathbf{A}$ is a known matrix of coefficients, $\mathbf{b}$ is a known right-hand side vector, and $\mathbf{x}$ is the unknown vector we want to find.

If $\mathbf{A}$ is square and invertible, the solution is simply:

However, there is one important practical point: never solve a linear system by explicitly computing the inverse. Computing $\mathbf{A}^{-1}$ and then multiplying is slower and numerically less stable than directly solving the system. Use a solver instead.

import numpy as np


# 2x +  y = 8
#  x + 3y = 9

A = np.array([[2, 1],
              [1, 3]])
b = np.array([8, 9])

# Preferred: solve directly
x = np.linalg.solve(A, b)
print(f"Solution: x = {x[0]}, y = {x[1]}")

# Check by multiplying
print(np.allclose(A @ x, b))  # True

# Compare with inverse (do not use in practice)
A_inv = np.linalg.inv(A)
x_inv = A_inv @ b
print(f"Same result via inverse: {x_inv}")

PROMPT

Use numpy to solve a 2x2 linear system with np.linalg.solve and verify the solution with matrix-vector multiplication.

A Financial Example: Portfolio Weights from Target Returns

Suppose you have a covariance matrix $\mathbf{\Sigma}$ and you want to find the portfolio weights $\mathbf{w}$ that produce a target expected return $\mu^*$ given that the expected return of each asset is $\boldsymbol{\mu}$. This leads to a system of linear equations:

With the constraint that weights sum to 1 ($\mathbf{w}^\top \mathbf{1} = 1$), you have two equations and $n$ unknowns. For a specific case with two assets, you can solve this directly:

import numpy as np

# Two assets
mu = np.array([0.08, 0.12])   # expected returns
ones = np.array([1.0, 1.0])    # sum-to-1 constraint

# Target portfolio return of 10%
target_return = 0.10

# Build the system A @ w = b
# Equation 1: w1 + w2 = 1
# Equation 2: 0.08*w1 + 0.12*w2 = 0.10
A = np.array([[1, 1],
              [mu[0], mu[1]]])
b = np.array([1, target_return])

w = np.linalg.solve(A, b)
print(f"Weights: {w}")
print(f"Sum: {w.sum():.4f}, Return: {w @ mu:.4f}")

PROMPT

Use numpy to solve a small linear system that finds portfolio weights summing to 1 and hitting a target return.

For larger systems with many assets and constraints, the same principles apply, you are solving $\mathbf{A}\mathbf{x} = \mathbf{b}$, and numpy.linalg.solve (or the QR decomposition behind it) is the right tool.

References

Strang, Introduction to Linear Algebra, 5th edition, Wellesley-Cambridge Press, Chapters 2–3. ISBN: 978-0980232776.

Lax, Linear Algebra and Its Applications, 2nd edition, Wiley, Chapters 1–2. ISBN: 978-0471751564.