Story Opening

The profiler’s verdict was unambiguous. The hot spot was this function, called once per transaction:

import statistics


def zscore_loop(amount: float, history: list[float]) -> float:
    mean = statistics.fmean(history)
    std = statistics.pstdev(history)
    return (amount - mean) / std if std else 0.0


print(round(zscore_loop(900.0, [100.0, 120.0, 80.0, 110.0, 90.0]), 2))  # -> 56.57

Called a hundred million times, it took most of an hour. In Java, the JIT would have compiled this into tight machine code and Arjun would have shrugged. In CPython, every iteration of every loop goes through the bytecode interpreter: each + looks up types, unboxes objects, allocates a new float object, and decrements reference counts.

Priya rewrote it with NumPy. The whole batch — every transaction against every customer’s history — ran in under a second.

“Python is the steering wheel,” she said. “NumPy is the engine. Your job is to make sure the loops happen in C, not in Python.”

NumPy is the foundation of the scientific Python stack. pandas is built on it, scikit-learn takes and returns its arrays, and PyTorch tensors are deliberately NumPy-shaped. Learn to think in arrays here, and every later library feels familiar.

Java → Python: The Quick Map

JavaNumPy
double[], double[][]np.ndarray (any number of dimensions)
new double[n]np.zeros(n) / np.empty(n)
for (i...) c[i] = a[i] + b[i]c = a + b
Arrays.stream(a).sum()a.sum()
Arrays.sort(a)np.sort(a) / a.sort()
System.arraycopy / subList viewSlicing returns a view
java.util.Random(seed)np.random.default_rng(seed)
Primitive overflow wraps silentlyFixed-width dtypes overflow too (int8, int32)
Vector API (incubator), EJML, ND4JNumPy — the standard everyone builds on

Deep Dive: Why NumPy Is Fast

A Python list of floats is an array of pointers, each pointing to a separate heap-allocated float object (24 bytes each, scattered in memory). A NumPy array is a single contiguous block of raw machine values with one shared type — exactly like a Java double[].

graph TB subgraph "Python list [1.0, 2.0, 3.0]" L[list header] --> P1[ptr] & P2[ptr] & P3[ptr] P1 --> F1["float obj 1.0"] P2 --> F2["float obj 2.0"] P3 --> F3["float obj 3.0"] end subgraph "np.array([1.0, 2.0, 3.0])" H["ndarray header
dtype=float64, shape=(3,)"] --> B["1.0 | 2.0 | 3.0
(24 contiguous bytes)"] end

Operations on arrays (a + b, a.sum(), np.sqrt(a)) run as compiled C loops over that buffer — cache-friendly, SIMD-vectorised, with no per-element interpreter overhead. This is called vectorisation: express the operation on the whole array and let NumPy run the loop.

import time
import numpy as np


n = 1_000_000
values = list(range(n))
arr = np.arange(n, dtype=np.float64)


start = time.perf_counter()
squared_list = [v * v for v in values]          # one million interpreter iterations
loop_s = time.perf_counter() - start


start = time.perf_counter()
squared_arr = arr * arr                         # one C loop
vec_s = time.perf_counter() - start


print(f"list comprehension: {loop_s * 1000:.1f} ms")
print(f"numpy vectorised:   {vec_s * 1000:.1f} ms")
print(f"speed-up: ~{loop_s / vec_s:.0f}x")      # typically 10-100x on a laptop

Tip — the golden rule: if you see a Python for loop iterating over the elements of an array, there’s almost always a vectorised way to write it. Loops over a handful of columns, files or epochs are fine.

Creating Arrays

import numpy as np


a = np.array([120.0, 45.5, 3000.0])          # from a list
print(a.dtype, a.shape, a.ndim)              # -> float64 (3,) 1


print(np.zeros(3))                           # -> [0. 0. 0.]
print(np.ones((2, 3), dtype=np.int32))       # 2 rows x 3 cols of int32 ones
print(np.full(3, 7.5))                       # -> [7.5 7.5 7.5]
print(np.arange(0, 10, 3))                   # -> [0 3 6 9]       (like range, but an array)
print(np.linspace(0, 1, 5))                  # -> [0.   0.25 0.5  0.75 1.  ]  (n evenly spaced points)
print(np.eye(2))                             # 2x2 identity matrix


# Reshape without copying data: 6 elements -> 2 rows x 3 columns
m = np.arange(6).reshape(2, 3)
print(m)
# [[0 1 2]
#  [3 4 5]]
print(m.reshape(-1, 2).shape)                # -> (3, 2)   (-1 = "infer this dimension")
print(m.T.shape)                             # -> (3, 2)   (transpose)
print(m.ravel())                             # -> [0 1 2 3 4 5]  (flatten)

dtypes — fixed-width types are back

NumPy arrays have one element type. The defaults are float64 and int64, but ML code often uses float32 to halve memory and match GPU expectations.

import numpy as np


x = np.array([1, 2, 3])
print(x.dtype)                               # -> int64
print((x / 2).dtype)                         # -> float64   (true division promotes)
print(x.astype(np.float32).nbytes)           # -> 12        (3 values x 4 bytes)


# GOTCHA: fixed-width integers OVERFLOW — just like Java's int, and unlike Python's int.
small = np.array([120, 127], dtype=np.int8)
print(small + 10)                            # -> [-126 -119]   (wrapped around!)


# Mixed input gets upcast to a common type.
print(np.array([1, 2.5]).dtype)              # -> float64
print(np.array([1, "a"]).dtype)              # -> <U21  (strings — almost always a bug)

Gotcha — If a CSV column contains a stray string, NumPy (and pandas) may give you an array of strings or Python objects (dtype=object), and every operation becomes slow or fails. Check .dtype early and often.

Vectorised Operations and Universal Functions

import numpy as np


amounts = np.array([120.0, 45.5, 3000.0, 80.0])
fx = np.array([1.0, 1.0, 83.0, 90.0])        # per-transaction FX rate


print((amounts * fx).tolist())               # -> [120.0, 45.5, 249000.0, 7200.0]  (element-wise)
print(amounts + 10)                          # scalar applied to every element
print(np.log1p(amounts).round(2))            # -> [4.8  3.84 8.01 4.39]  (ufunc: log(1 + x))
print(np.sqrt(np.array([4.0, 9.0])))         # -> [2. 3.]
print(np.clip(amounts, 50, 1000))            # -> [ 120.   50. 1000.   80.]


# Aggregations
print(amounts.sum(), amounts.mean(), amounts.max())   # -> 3245.5 811.375 3000.0
print(amounts.std().round(2))                         # -> 1263.88   (population std, ddof=0)
print(amounts.argmax())                               # -> 2   (index of the max)
print(np.percentile(amounts, 50))                     # -> 100.0   (median)
print(np.cumsum([1, 2, 3]))                           # -> [1 3 6]

Gotcha — std defaults differ. NumPy’s std() is the population standard deviation (ddof=0). pandas’ .std() is the sample standard deviation (ddof=1). Same data, different answers. Pass ddof explicitly when it matters.

Indexing: Slices, Masks and Fancy Indexing

import numpy as np


# Rows = customers, columns = daily spend over 5 days
spend = np.array([
    [100, 120,  80, 110,  90],
    [ 20,  25,  30,  20, 900],
    [500, 520, 480, 510, 495],
])


print(spend[1, 4])          # -> 900          (row 1, col 4 — one bracket, comma-separated)
print(spend[0])             # -> [100 120  80 110  90]   (first row)
print(spend[:, -1])         # -> [ 90 900 495]          (last column, all rows)
print(spend[:2, 1:3])       # rows 0-1, cols 1-2 -> [[120 80] [25 30]]


# Boolean masks: a comparison produces an array of booleans...
big = spend > 400
print(big[1])               # -> [False False False False  True]
# ...which you can use to SELECT elements (result is 1-D):
print(spend[big])           # -> [900 500 520 480 510 495]
print(big.sum())            # -> 6    (True counts as 1: "how many?")


# Combine conditions with & | ~ and PARENTHESES (not 'and'/'or'/'not').
print(spend[(spend > 100) & (spend < 500)])   # -> [120 110 480 495]


# np.where: vectorised ternary
print(np.where(spend[1] > 100, "spike", "ok"))  # -> ['ok' 'ok' 'ok' 'ok' 'spike']


# Fancy indexing: select with an array of indices
print(spend[[2, 0], 0])     # -> [500 100]   (rows 2 and 0, column 0)

Gotcha — and / or don’t work on arrays. (a > 1) and (a < 5) raises ValueError: The truth value of an array with more than one element is ambiguous. Python’s and needs a single bool; use the element-wise operators &, |, ~ — and always parenthesise each comparison, because & binds tighter than >.

Deep Dive: Axes

Aggregations take an axis argument, and it confuses everyone at first. The trick: axis is the dimension that gets collapsed (summed away).

For a 2-D array of shape (rows, cols):

  • axis=0 collapses the rows → one result per column (shape (cols,)).
  • axis=1 collapses the columns → one result per row (shape (rows,)).
import numpy as np


spend = np.array([
    [100, 120,  80],     # customer A
    [ 20,  25,  30],     # customer B
])                       # shape (2 customers, 3 days)


print(spend.sum())               # -> 375          (everything)
print(spend.sum(axis=0))         # -> [120 145 110]  (per DAY: rows collapsed)
print(spend.sum(axis=1))         # -> [300  75]     (per CUSTOMER: columns collapsed)
print(spend.mean(axis=1, keepdims=True))   # keepdims keeps a (2, 1) shape — useful for broadcasting
# [[100.]
#  [ 25.]]

In ML, the convention is almost universal: rows are samples, columns are features. So “normalise each feature” means aggregating with axis=0, and “score each sample” means axis=1.

Deep Dive: Broadcasting

Broadcasting is how NumPy combines arrays of different shapes without copying data. You’ve already used it: amounts + 10 broadcasts the scalar 10 across the array. The general rule:

Compare shapes from the right. Two dimensions are compatible if they are equal or one of them is 1. A missing dimension counts as 1. Size-1 dimensions are stretched (virtually) to match.

Shape AShape BResultWhy
(3,)() scalar(3,)scalar stretches everywhere
(4, 3)(3,)(4, 3)B treated as (1, 3), stretched over 4 rows
(4, 3)(4, 1)(4, 3)column vector stretched across 3 columns
(4, 1)(1, 3)(4, 3)both stretch — an “outer” operation
(4, 3)(4,)error(4,) aligns with the 3 on the right: 3 ≠ 4

Here is Arjun’s z-score problem solved for every customer at once. Each customer’s latest transaction is compared against their own history — no loops:

import numpy as np


# history: 3 customers x 5 past transactions
history = np.array([
    [100.0, 120.0,  80.0, 110.0,  90.0],
    [ 20.0,  25.0,  30.0,  20.0,  25.0],
    [500.0, 520.0, 480.0, 510.0, 490.0],
])
latest = np.array([900.0, 26.0, 505.0])     # one new transaction per customer, shape (3,)


mean = history.mean(axis=1)                 # shape (3,)  — per customer
std = history.std(axis=1)                   # shape (3,)


z = (latest - mean) / std                   # (3,) op (3,) op (3,) -> (3,)
print(z.round(2))                           # -> [56.57  0.53  0.35]


# Now z-score EVERY historical value against its own customer's stats.
# history is (3, 5); mean is (3,). Shapes compared from the right: 5 vs 3 -> ERROR.
try:
    history - mean
except ValueError as e:
    print("broadcast error:", "could not be broadcast" in str(e))  # -> broadcast error: True


# Fix: make mean a COLUMN vector of shape (3, 1) so it stretches across the 5 columns.
z_all = (history - mean[:, np.newaxis]) / std[:, np.newaxis]
# Equivalent: history.mean(axis=1, keepdims=True)
print(z_all.shape)                          # -> (3, 5)
print(np.allclose(z_all.mean(axis=1), 0))   # -> True   (each row now has mean 0)

Feature standardisation for ML — “subtract each column’s mean, divide by its std” — is the same idea along the other axis, and it needs no reshaping because (n, features) - (features,) broadcasts naturally:

import numpy as np


rng = np.random.default_rng(seed=42)
X = rng.normal(loc=[100, 3, 0.5], scale=[30, 1, 0.1], size=(1_000, 3))   # 1000 samples, 3 features


X_std = (X - X.mean(axis=0)) / X.std(axis=0)    # (1000, 3) - (3,) -> broadcasts per column
print(X_std.mean(axis=0).round(6) + 0.0)        # -> [0. 0. 0.]
print(X_std.std(axis=0).round(6))               # -> [1. 1. 1.]

That two-line function is what scikit-learn’s StandardScaler does internally (Part 9).

Tip — when shapes confuse you, print them. print(a.shape, b.shape) before an operation resolves 90% of broadcasting bugs. The most dangerous case is the one that doesn’t error: (n,) vs (n, 1) silently broadcasts to (n, n) — a giant matrix of wrong answers.

import numpy as np


y_true = np.array([1.0, 0.0, 1.0])          # shape (3,)
y_pred = np.array([[0.9], [0.2], [0.8]])    # shape (3, 1) — e.g. a model output


print((y_true - y_pred).shape)              # -> (3, 3)   (SILENT BUG: not element-wise!)
print((y_true - y_pred.ravel()).shape)      # -> (3,)     (correct)

Deep Dive: Views vs Copies

Basic slicing returns a view: a new array header pointing into the same memory. It’s fast (no copying) — and it means writes through the view change the original. Think of Java’s List.subList(), but for every slice.

import numpy as np


prices = np.array([10.0, 20.0, 30.0, 40.0])


window = prices[1:3]                 # VIEW — shares memory
window[0] = 999.0
print(prices)                        # -> [ 10. 999.  30.  40.]   (original changed!)
print(np.shares_memory(prices, window))   # -> True


# Boolean and fancy indexing return COPIES.
selected = prices[prices > 25]
selected[:] = 0
print(prices)                        # -> [ 10. 999.  30.  40.]   (unchanged)


# Need an independent slice? Copy explicitly.
safe = prices[1:3].copy()
print(np.shares_memory(prices, safe))     # -> False


# Assigning THROUGH a mask on the original modifies it in place (common and useful):
prices[prices > 100] = 100.0
print(prices)                        # -> [ 10. 100.  30.  40.]

Gotcha — This is the opposite of Python lists, where slices are copies (Part 2). It’s deliberate: copying gigabytes on every slice would make NumPy unusable. Remember: slices are views; masks and index arrays are copies. reshape and .T also return views when possible.

Missing Values: NaN

NumPy represents missing floats as np.nan (IEEE “not a number”). NaN poisons arithmetic and isn’t equal to anything — including itself.

import numpy as np


amounts = np.array([120.0, np.nan, 80.0])


print(amounts.mean())               # -> nan       (NaN propagates)
print(np.nanmean(amounts))          # -> 100.0     (nan-aware variants: nanmean, nansum, nanmax...)
print(np.nan == np.nan)             # -> False     (never compare to NaN with ==)
print(np.isnan(amounts))            # -> [False  True False]
print(np.nan_to_num(amounts, nan=0.0))   # -> [120.   0.  80.]

pandas builds a richer missing-data story on top of this (Part 8).

Random Numbers, Reproducibly

import numpy as np


# The modern API: a Generator object with an explicit seed (≈ new Random(42)).
rng = np.random.default_rng(seed=42)


print(rng.integers(0, 10, size=5))            # 5 ints in [0, 10)
amounts = rng.lognormal(mean=4.0, sigma=1.0, size=10_000)   # skewed, like real payments
print(amounts.min() > 0)                      # -> True
is_fraud = rng.random(10_000) < 0.02          # ~2% positives: a boolean mask
print(0.01 < is_fraud.mean() < 0.03)          # -> True


idx = rng.permutation(10)                     # shuffled indices — for manual train/test splits
print(sorted(idx) == list(range(10)))         # -> True


# Same seed, same numbers — essential for reproducible experiments.
print((np.random.default_rng(7).random(3) == np.random.default_rng(7).random(3)).all())  # -> True

Tip — Avoid the legacy global API (np.random.seed(...), np.random.rand(...)) in new code. A passed-around Generator makes randomness explicit and keeps components from interfering with each other’s streams.

Linear Algebra: Similarity Search From Scratch

Embeddings — vectors that represent text, images or customers — are central to modern AI. Finding “similar” items means comparing vectors, usually by cosine similarity: the dot product of two vectors divided by the product of their lengths. Here is a complete, vectorised nearest-neighbour search over 10,000 merchant embeddings:

import numpy as np


rng = np.random.default_rng(seed=0)


n_merchants, dim = 10_000, 64
embeddings = rng.normal(size=(n_merchants, dim)).astype(np.float32)   # (10000, 64)


# Make merchant 7 a near-duplicate of merchant 42 so we know the right answer.
embeddings[7] = embeddings[42] + rng.normal(scale=0.05, size=dim)


query = embeddings[42]                                                 # (64,)


# 1. Normalise every vector to unit length (axis=1: per row), keepdims for broadcasting.
norms = np.linalg.norm(embeddings, axis=1, keepdims=True)              # (10000, 1)
unit = embeddings / norms                                              # (10000, 64)
q = query / np.linalg.norm(query)                                      # (64,)


# 2. Cosine similarity of the query against ALL merchants: one matrix-vector product.
scores = unit @ q                                                      # (10000,)  '@' = matmul


# 3. Top-k without fully sorting: argpartition is O(n); then sort just those k.
k = 3
top = np.argpartition(-scores, k)[:k]
top = top[np.argsort(-scores[top])]


print(top[:2])                                    # -> [42  7]
print(np.round(scores[top[0]], 4))                # -> 1.0

That’s the core of a vector database, in about ten lines. Production systems (FAISS, pgvector, Redis vector search) add approximate indexes like HNSW so they don’t have to scan every row.

Operator / functionMeaning
a * bElement-wise product
a @ b / np.matmulMatrix multiplication (dot product for 1-D)
a.TTranspose
np.linalg.norm(a, axis=...)Vector length(s)
np.linalg.inv, solve, eig, svdInverse, linear systems, eigen-decomposition, SVD

Tips, Tricks & Gotchas

Gotcha — growing arrays in a loop. arr = np.append(arr, x) copies the whole array every time: O(n²). Collect values in a Python list and convert once with np.array(values), or preallocate with np.empty(n) and fill by index.

Tip — np.vectorize is not vectorisation. It’s a convenience wrapper around a Python loop with the same speed. Look for a real ufunc, np.where, or array arithmetic instead.

Tip — readable printing. np.set_printoptions(precision=3, suppress=True) stops scientific notation from cluttering your output in notebooks.

Gotcha — integer division of int arrays. np.array([7]) // 2 is floor division; / returns floats. Same rules as plain Python (Part 1), applied element-wise.

Tip — memory math. A float64 matrix of 10 million rows × 50 features is 4 GB. Switching to float32 halves it. Always do the arithmetic before loading “all the data” into memory.

Key Takeaways

ConceptRemember
ndarrayContiguous, typed, n-dimensional; operations run in C
VectorisationReplace element loops with whole-array expressions
dtypesfloat64 default, float32 for ML; fixed-width ints overflow
Axesaxis = the dimension collapsed; rows = samples, columns = features
Masksa[(a > x) & (a < y)]; use &, |, ~ with parentheses
BroadcastingAlign shapes from the right; equal or 1; watch (n,) vs (n, 1)
Views vs copiesSlices are views; masks and fancy indexes are copies
NaNPropagates; use nan* functions and np.isnan
Randomnessrng = np.random.default_rng(seed)

Story Closing

The vectorised z-score ran over a hundred million transactions in 1.8 seconds. Arjun ran it three times because he didn’t believe it.

Flush with success, he tried to build the rest of the feature pipeline in raw NumPy: per-merchant averages, time-of-day buckets, joins against the customer table. Within an hour he was juggling seven parallel arrays, keeping their indices aligned by hand, and storing column meanings in comments. One off-by-one in a sort and every feature would silently belong to the wrong customer.

“You’ve reinvented a table,” Priya said. “Badly. You want labelled columns, group-bys and joins. You want pandas — think of it as SQL inside Python.”

In Part 8, Arjun learns pandas through the lens he already knows best: SQL.

This is Part 7 of a 10-part series: “Python for Java Developers: From Streams to Tensors.”