Vectorisation

Think in whole arrays, not in elements

Vectorisation means writing a calculation as operations on whole arrays, and letting NumPy run the loop in fast compiled code. a + b adds two arrays element by element, a * b multiplies them, a ** 2 squares every element and np.sqrt(a) takes every square root. These are universal functions (ufuncs). The result is shorter to write, easier to read, and, for large data, often fifty times faster than the equivalent Python for loop, which pays interpreter overhead on every single element.

import numpy as np

a = np.array([1.0, 2.0, 3.0])
b = np.array([10.0, 20.0, 30.0])

print(a + b)             # [11. 22. 33.]
print(a * b)             # [10. 40. 90.]
print(a ** 2)            # [1. 4. 9.]
print(np.sqrt(a))        # [1.         1.41421356 1.73205081]

# The loop version of a * b, which vectorisation replaces
result = np.empty(3)
for i in range(3):
    result[i] = a[i] * b[i]
print(result)            # [10. 40. 90.]

Returns without a loop: shifting with slices

A calculation that compares each element with its neighbour needs a second copy of the array shifted by one, and slices give it for free. p[1:] is today's price and p[:-1] yesterday's, both of length n - 1, so the simple returns are p[1:] / p[:-1] - 1. The same pair gives differences (p[1:] - p[:-1], or np.diff(p)). Running totals and products are np.cumsum and np.cumprod; a growth path is np.cumprod(1 + returns). Whole-array comparisons (r > 0) and np.where(condition, if_true, if_false) replace if inside a loop.

import numpy as np

p = np.array([100.0, 102.0, 101.0, 103.0])

returns = p[1:] / p[:-1] - 1                 # today / yesterday - 1, no loop
print(np.round(returns, 4))                  # [ 0.02   -0.0098  0.0198]
print(np.diff(p))                            # [ 2. -1.  2.]
print(np.cumprod(1 + np.array([0.01, 0.02])))   # [1.01   1.0302]

# np.where replaces if/else inside a loop
r = np.array([0.01, -0.02, 0.03, -0.01])
print(np.where(r > 0, "up", "down"))         # ['up' 'down' 'up' 'down']
print(np.where(r > 0, r, 0.0))               # [0.01 0.   0.03 0.  ]

Aggregations along an axis, and matrix products

Functions such as sum, mean, std, min and max collapse an array. With a 2-D array of returns (one row per day, one column per asset), axis=0 collapses the rows and gives one value per column (each asset's average over time), while axis=1 gives one value per row (each day's average across assets). The matrix product R @ w, of a (days, assets) array with a weight vector of length assets, gives the portfolio return of every day in one call: that is the dot product from the vectors topic done for all rows at once.

import numpy as np

# 2 days x 3 assets
R = np.array([[0.02, -0.01, 0.03],
              [0.01,  0.00, 0.02]])
w = np.array([0.5, 0.3, 0.2])

print(R.shape)            # (2, 3)
print(R.mean(axis=0))     # [ 0.015 -0.005  0.025]   one value per asset
print(R.sum(axis=1))      # [0.04 0.03]              one value per day
print(R @ w)              # [0.013 0.009]            portfolio return each day

# Standardise: how many standard deviations from the mean?
z = np.array([1.0, 2.0, 3.0, 4.0])
print((z - z.mean()) / z.std())   # [-1.34164079 -0.4472136   0.4472136   1.34164079]

Moving averages with cumsum, and the limits of vectorising

A moving average over a window w looks like it needs a loop, but with a running total it does not: put a 0 in front of np.cumsum(x) and the sum of any w consecutive values is a difference of two running totals, cs[w:] - cs[:-w]; divide by w. The result has len(x) - w + 1 values. Vectorise the arithmetic, not logic that depends on the previous result (a position that depends on yesterday's fill, say): those stay in a loop, or move to compiled tools such as Numba later. Watch memory too, because every intermediate array is a full copy.

import numpy as np

x = np.array([1.0, 2.0, 4.0, 8.0, 16.0])
w = 3

cs = np.concatenate(([0.0], np.cumsum(x)))   # [0, 1, 3, 7, 15, 31]
print((cs[w:] - cs[:-w]) / w)                # [2.33333333 4.66666667 9.33333333]

# Check against the plain definition
print([x[i:i + w].mean() for i in range(len(x) - w + 1)])   # three values, the same

• Whole-array operations (a + b, a ** 2, np.sqrt(a)) replace loops; compare with np.where.
• Shifted slices: p[1:] / p[:-1] - 1; cumsum and cumprod for running totals and growth.
• axis=0 collapses rows (per column), axis=1 collapses columns (per row); R @ w is a dot product per row.

Next up: broadcasting, the rule that lets arrays of different shapes work together.

Lab · what's graded

• Compute all the returns at once with simple_returns25 pts
• Get the portfolio return of every day with portfolio_returns25 pts
• Smooth a series without a loop using moving_average25 pts
• Standardise values with zscore25 pts