NumPy & pandas

Every numerical library in this field expects whole-array operations rather than loops, and the gap is not small: the loop version of a job can be a hundred times slower than the same thing expressed as one expression over the whole array. This section is about thinking that way, along with the indexing, masking and reshaping that make it possible, and the table operations in pandas that most real data work actually consists of.

Stop writing loops#

A loop over a NumPy array is usually a hundred times slower than the vectorised equivalent, because the loop runs in Python while the vectorised version runs in compiled code. Learning to see operations as whole-array is the main skill here.

import numpy as np

# Slow, and how most people first write it
result = []
for x in data:
    result.append(x * 2 + 1)

# Fast, and shorter
result = data * 2 + 1

Two things make this worth the adjustment. The obvious one is speed: the Python version pays interpreter overhead per element, boxing each number into an object, while the array version dispatches once into compiled code that walks contiguous memory. The less obvious one is that every library you will hand this data to — scikit-learn, PyTorch, anything — expects arrays, so a list of Python floats has to be converted anyway.

The axis argument#

This is the single thing people get wrong most often, and the rule is short: axis names the dimension that disappears. An array of shape (2, 3) summed along axis 0 loses the 2 and leaves shape (3,); summed along axis 1 it loses the 3 and leaves shape (2,). Say it that way rather than in terms of rows and columns, because the rows-and-columns version stops working the moment there are three dimensions.

a = np.array([[1, 2, 3],
              [4, 5, 6]])            # shape (2, 3)

a.sum()                              # 21          everything, shape ()
a.sum(axis=0)                        # [5 7 9]     shape (3,)  the 2 is gone
a.sum(axis=1)                        # [ 6 15]     shape (2,)  the 3 is gone
a.sum(axis=1, keepdims=True).shape   # (2, 1)      kept as a length-1 dimension
a.argmax(axis=1)                     # [2 2]       index of the max in each row

The same rule applies to any array. For a batch of token embeddings of shape (batch, sequence, features), averaging over axis=1 collapses the sequence and gives one vector per item, which is mean pooling. Averaging over axis=0 collapses the batch, which is almost never what anyone meant. Negative indices count from the end, so axis=-1 is the last dimension and is the safe way to write an operation that should apply to features regardless of how many leading dimensions there are.

def softmax(z):
    """Over the last axis, so it works on (n,), (batch, n) and (batch, seq, n)."""
    z = z - z.max(axis=-1, keepdims=True)      # subtract the max: stops exp overflowing
    e = np.exp(z)
    return e / e.sum(axis=-1, keepdims=True)

logits = np.array([[2.0, 1.0, 0.1],
                   [0.5, 0.5, 0.5]])
softmax(logits)
# [[0.659 0.242 0.099]
#  [0.333 0.333 0.333]]
softmax(logits).sum(axis=-1)         # [1. 1.]  always check this

Broadcasting, as rules#

Broadcasting is how NumPy lets arrays of different shapes take part in the same operation without you writing a loop or copying data. It is worth learning as three mechanical rules rather than as intuition, because then you can work out what a line does without running it.

1. Line the shapes up from the RIGHT, padding the shorter one with 1s on the left.
2. Two dimensions are compatible if they are equal, or if one of them is 1.
3. A dimension of 1 is stretched to match the other. Anything else is an error.

  (3, 2)      (3, 1)        (2, 3)        (3, 2)
+ (   2)    + (   2)      + (2, 1)      + (   3)
---------   ---------     --------      ---------
  (3, 2)      (3, 2)        (2, 3)        error: 2 and 3 are neither
   ok          ok            ok           equal nor 1

x = np.arange(6).reshape(3, 2)       # [[0 1] [2 3] [4 5]]
row = np.array([10, 20])             # (2,)   -> treated as (1, 2), stretched down rows
col = np.array([[100], [200], [300]])  # (3, 1) -> stretched across columns

x + row
# [[10 21]
#  [12 23]
#  [14 25]]

x + col
# [[100 101]
#  [202 203]
#  [304 305]]

x + np.array([1, 2, 3])
# ValueError: operands could not be broadcast together with shapes (3,2) (3,)

Nothing is copied. A stretched dimension is read repeatedly from the same memory, which is why subtracting a mean vector from a million rows costs nothing extra in space. The tool for making shapes line up is None, also spelled np.newaxis, which inserts a dimension of length 1 wherever you put it.

v = np.array([1, 2, 3])
v[:, None].shape      # (3, 1)  a column
v[None, :].shape      # (1, 3)  a row

# All pairwise differences, with no loop:
v[:, None] - v[None, :]
# [[ 0 -1 -2]
#  [ 1  0 -1]
#  [ 2  1  0]]

# Centre every feature of a (n, d) matrix on its own mean:
X - X.mean(axis=0)                 # (n,d) - (d,)  -> fine, stretches down rows

# Scale every row to unit length:
X / np.linalg.norm(X, axis=1, keepdims=True)    # (n,d) / (n,1)

Masking#

A boolean array selects elements. This replaces most filtering loops and composes cleanly.

scores = np.array([0.2, 0.9, 0.5, 0.95])

scores > 0.8                  # array([False, True, False, True])
scores[scores > 0.8]          # array([0.9, 0.95])
np.where(scores > 0.8, 1, 0)  # array([0, 1, 0, 1])

(scores > 0.4) & (scores < 0.9)   # & and |, not and/or

Three things follow from this. A mask is itself an array, so masks combine with & and | and invert with ~, and each operand needs its own brackets because & binds tighter than the comparison operators. A mask counts as well as selects, because summing booleans counts the Trues. And indexing with a mask, or with a list of indices, returns a copy rather than a view, which is the opposite of what a plain slice does.

mask = (scores > 0.4) & (scores < 0.9)

mask.sum()                # 1     how many match
mask.mean()               # 0.25  what fraction match
scores[~mask]             # everything that does not
np.where(mask)[0]         # array([2])  the positions, not the values

np.shares_memory(scores, scores[1:3])      # True:  a slice is a view
np.shares_memory(scores, scores[mask])     # False: a mask gives a copy
np.shares_memory(scores, scores[[0, 2]])   # False: so does a list of indices

scores[scores > 0.8] = 0.8    # assignment through a mask does write back

Shapes#

a = np.arange(12)

a.reshape(3, 4)         # explicit
a.reshape(3, -1)        # -1 means "work it out"
a[:, None].shape        # (12, 1)  add an axis
a.T                     # transpose
np.concatenate([a, a])  # join along an existing axis
np.stack([a, a])        # join along a new axis

reshape never moves data; it reinterprets the same buffer, so it is free and it returns a view. It fills the new shape in row-major order, last axis fastest, which is why reshaping (2, 3) to (3, 2) is not a transpose and gives a different arrangement of the same numbers. Transposing also returns a view, one whose strides are swapped rather than whose memory is rearranged, which is why some later operation on it may quietly copy.

a = np.arange(6).reshape(2, 3)
# [[0 1 2]
#  [3 4 5]]

a.reshape(3, 2)      # [[0 1] [2 3] [4 5]]   filled across, then down
a.T                  # [[0 3] [1 4] [2 5]]   not the same thing

np.concatenate([a, a], axis=0).shape   # (4, 3)  more rows
np.concatenate([a, a], axis=1).shape   # (2, 6)  more columns
np.stack([a, a]).shape                 # (2, 2, 3)  a new leading axis

a.squeeze()          # drop every axis of length 1
a.ravel()            # flatten to 1-D, a view when it can be

pandas#

A DataFrame is a table with labelled columns. The operations that matter are selection, grouping and joining, which cover most data preparation.

import pandas as pd

df = pd.read_csv("orders.csv")

df.head()                       # look at it, always
df.info()                       # types and null counts
df.describe()                   # distributions

df[df.total > 100]              # filter rows
df[["id", "total"]]             # select columns
df.loc[df.total > 100, "flag"] = True   # conditional assignment

df.groupby("customer_id").agg(
    orders=("id", "count"),
    revenue=("total", "sum"),
)

pd.merge(orders, customers, on="customer_id", how="left")

Two pieces of vocabulary before anything else. A DataFrame has an index as well as columns, and the index is what alignment happens on: adding two Series lines them up by index label, not by position, so a subtraction between a filtered frame and an unfiltered one produces NaN wherever the labels do not match rather than the answer you expected. And there are two accessors: .loc selects by label and .iloc by position, and mixing them up on a frame whose index is not 0..n-1 gives wrong rows silently.

df.loc[3, "total"]          # the row whose index LABEL is 3
df.iloc[3]                  # the FOURTH row, whatever its label
df.loc[df.total > 100, "flag"] = True     # label-based, and the way to assign

sub = df[df.total > 100]    # the index keeps the original labels: 4, 9, 17...
sub.reset_index(drop=True)  # renumber when position is what you want

The related habit is to call .copy() when you deliberately want an independent frame to work on, as in ok = df[df.status == 'ok'].copy(). It states the intention, and it means later columns you add are yours rather than an attempt to write through a filter.

Group-by and merge, the two that matter#

Group-by splits the frame by a key, applies something to each piece, and combines the results. Which of the three shapes you want depends on what comes out: agg gives one row per group, transform gives one value per original row, and apply gives whatever your function returns and is the slow one to avoid when either of the others will do.

calls = pd.DataFrame({
    "user_id":    [1, 1, 2, 2, 3],
    "model":      ["sonnet", "haiku", "sonnet", "sonnet", "haiku"],
    "tokens":     [1200, 300, 900, 1500, None],
    "latency_ms": [820, 140, 760, 2100, 155],
})

calls.groupby("model").agg(
    calls=("user_id", "size"),           # rows in the group
    users=("user_id", "nunique"),
    tokens=("tokens", "sum"),
    p95=("latency_ms", lambda s: s.quantile(0.95)),
)
#         calls  users  tokens      p95
# model
# haiku       2      2   300.0   154.25
# sonnet      3      2  3600.0  1972.00

Note haiku's tokens: two rows, one of them missing, and the sum is 300 rather than an error or a NaN. size counts rows in the group and count counts non-null values in a column, so on that group they give 2 and 1. Choosing the wrong one is how a denominator quietly becomes too small.

calls.groupby("model").size()          # haiku 2, sonnet 3   rows
calls.groupby("model")["tokens"].count()   # haiku 1, sonnet 3   non-null

# transform: one value per original row, so it can go straight back on the frame
calls["model_total"] = calls.groupby("model")["tokens"].transform("sum")
#    user_id   model  tokens  latency_ms  model_total
# 0        1  sonnet  1200.0         820       3600.0
# 1        1   haiku   300.0         140        300.0
# ...
calls["share"] = calls["tokens"] / calls["model_total"]

merge is a SQL join with a different spelling, and the same things go wrong. Two arguments turn its failure modes into errors instead of wrong numbers: validate states the relationship you believe holds and raises if it does not, and indicator adds a column saying which side each row came from so you can count the rows that failed to match.

users = pd.DataFrame({"user_id": [1, 2, 4], "plan": ["pro", "free", "pro"]})

m = calls.merge(users, on="user_id", how="left",
                validate="many_to_one", indicator=True)
m["_merge"].value_counts()
# both          4
# left_only     1        <- user 3 has no row in users
# right_only    0

# If users had two rows for user_id 1, the same call raises rather than silently
# duplicating: MergeError: Merge keys are not unique in right dataset;
# not a many-to-one merge
#
# Without validate, 5 rows in become 7 rows out and every sum is now wrong.

The habit worth forming is to print len(df) before and after every merge. A row count that changed when you expected a lookup is the single most common cause of a number being wrong in a data pipeline, and it is free to check.

Missing values#

df.isna().sum()                 # how many per column
df.dropna(subset=["total"])     # drop rows missing something essential
df.total.fillna(df.total.median())   # impute

# Whether a value is missing is often informative in itself
df["total_was_missing"] = df.total.isna()

This is where array code goes quietly wrong. Slicing usually returns a view onto the same memory rather than a copy, so writing to what you thought was a separate array changes the original — no error, wrong numbers. This section covers views and copies, choosing data types so a table fits in memory at all, what to do when it does not, and how to recognise the point where pandas has stopped being the right tool for the size of data in front of you.

Views and copies#

An array is a block of memory plus a description of how to walk it. A view is a second description of the same block, which is why slicing costs nothing regardless of how large the array is, and why writing through a slice writes into the original. This is efficient and it is the source of confusing bugs, because nothing in the syntax tells you which you have.

a = np.arange(10)
b = a[:5]         # a view
b[0] = 99
a[0]              # 99  -- the original changed

c = a[:5].copy()  # an independent copy

The rule is short. Basic indexing — a slice, a step, a transpose, a reshape — returns a view. Advanced indexing — a boolean mask, a list or array of indices — returns a copy, because the selected elements are not evenly spaced and cannot be described by strides. Assignment through either does write to the original, which is the part that makes the rule feel inconsistent until you see why.

np.shares_memory(a, a[2:7])      # True   slice: a view
np.shares_memory(a, a[::2])      # True   step: still a view
np.shares_memory(a, a.T)         # True   transpose: a view
np.shares_memory(a, a[a > 3])    # False  mask: a copy
np.shares_memory(a, a[[0, 2]])   # False  list of indices: a copy

a[2:7].flags["OWNDATA"]          # False  it is borrowing
a[2:7].base is a                 # True   this is who it borrowed from

a[a > 3] = 0                     # but assignment through a mask DOES write back

pandas has the same question and has answered it differently over time. Under the copy-on-write behaviour that became the default in pandas 3, every operation behaves as though it returned a copy, and the copy is only actually made if you write to it. That removes the ambiguity: chained assignment now raises ChainedAssignmentError rather than the old SettingWithCopyWarning, and .loc is the way to assign in both.

Data types#

Memory scales directly with the data type, and the default is generous. A column of ages stored as the default 64-bit integer uses eight bytes per value to hold a number that fits in one. On a million rows that is 8 MB against 1 MB, and on a table with thirty such columns it is the difference between fitting in memory and not.

dtype      bytes   range or precision              use it for
int8         1      -128 to 127                     small counts, ages, flags
int16        2      -32,768 to 32,767               most counts
int32        4      about +/- 2.1 billion           ids, unless they exceed that
int64        8      the default                     large ids, timestamps in ns
float32      4      about 7 significant digits      embeddings, activations, images
float64      8      about 16 significant digits     the default; money, accumulators
bool         1      True or False
category   varies   stores each distinct string once  repeated strings in pandas

df.memory_usage(deep=True).sum()     # 70,500,000 on a 1M-row frame

df["age"] = df["age"].astype("int8")
df["country"] = df["country"].astype("category")

df.memory_usage(deep=True).sum()     # 2,000,000 -- 35x smaller

Categorical types for repeated strings are frequently the largest single saving, because a column of a million country names becomes a million small integers plus a lookup table of the handful of distinct values. Use it when the number of distinct values is small relative to the rows, and not when nearly every value is unique, where it adds overhead instead.

np.array([120, 100], dtype="int8") + np.int8(10)    # [-126  110], no warning
np.int8(120) + np.int8(10)                          # -126, with a RuntimeWarning

np.float32(1e8) + np.float32(1.0) == np.float32(1e8)    # True

arr = np.arange(5)                  # int64
arr += 0.5
# UFuncTypeError: Cannot cast ufunc 'add' output from dtype('float64')
#                 to dtype('int64') with casting rule 'same_kind'
arr = arr + 0.5                     # this works, and returns a new float64 array

Where the memory actually goes#

A MemoryError almost never happens at the line that allocated the data. It happens at some intermediate result: every arithmetic expression allocates a new array of the same size, so a chain of four operations on a 6 GB matrix briefly needs several copies of it. The fixes are in-place operators, which write into the existing buffer, and doing the work in chunks.

X = X - mean          # allocates a whole new X
X /= std              # in place: no allocation

np.subtract(X, mean, out=X)        # explicit form of the same thing

# Chunking, when even one pass will not fit:
total = 0.0
for chunk in pd.read_csv("events.csv", chunksize=500_000):
    total += chunk["cost_usd"].sum(dtype="float64")

Data larger than memory#

• Read in chunks and process each, so memory stays bounded. pandas read_csv takes chunksize and returns an iterator of frames.
• Use a columnar format such as Parquet and read only the columns you need. A twenty-column CSV read for two columns parses all twenty; Parquet reads two.
• Give read_csv the dtypes rather than letting it infer them. Inference reads the file to guess, then stores the widest type it saw, and a single stray string turns a numeric column into object.
• For larger analytical work, an engine like DuckDB or Polars handles data beyond memory efficiently and is often far faster than pandas on a single machine.
• The point at which pandas is the wrong tool is not a fixed size. It is when one machine's memory cannot hold a few copies of the data, which in practice is somewhere around a few gigabytes on a laptop and tens of gigabytes on a server.

pandas is where most data work in Python actually happens: a table held in memory, with named columns and typed values, that you filter, group, join and reshape. This section takes one real dataset from a CSV file through to a weekly summary per model, doing each step the way that survives contact with bad data. The habits are the point, particularly declaring types on load, looking at the rows before transforming them, and asking a join to check its own assumptions. They are what separates a number you can quote from one that is quietly wrong.

Step            Command                                    What it saves you from
load            read_csv(dtype=..., parse_dates=...)       silent object columns, 10x memory
look            info(), describe(), value_counts()         a whole analysis of bad data
filter          df[mask].copy()                            writing through a filtered view
derive          vectorised arithmetic, not apply           a job that takes an hour
aggregate       groupby().agg(named=("col", "fn"))         unreadable multi-level columns
join            merge(validate=..., indicator=True)        rows multiplying, totals doubling
resample        set_index(date).resample("D")              missing days vanishing from a chart

Load it properly#

import pandas as pd

calls = pd.read_csv(
    "llm_calls.csv",
    parse_dates=["created_at"],
    dtype={"model": "category", "status": "category", "user_id": "int64"},
)

Stating the types does three jobs. It stops the inference that turns a column with one stray value into object, which is pandas storing Python objects and is both slow and eight times the memory. It makes dates into real timestamps so that comparison, sorting and resampling work, rather than strings where '2026-3-1' sorts after '2026-12-01'. And a category column stores each distinct string once, which on a model or status column is usually the largest saving in the frame. Add usecols to read only the columns you need, and nrows=1000 while you are still working out what the file contains.

Look first#

calls.info()                           # columns, types, missing counts, memory
calls.describe()                       # numeric summaries: look for impossible values
calls["status"].value_counts(dropna=False)
calls.sample(5, random_state=0)        # read actual rows

Five minutes here saves hours later. Negative token counts, a status of 'OK ' with a trailing space, or a column that is 40% empty change everything you do next.

Read each of those for something specific. info tells you whether any column came in as object when you expected a number, and how many nulls each has. describe shows min and max, which is where impossible values live: a negative token count, a latency of zero, a date in 1970 that is really a missing value written as an epoch. value_counts with dropna=False is the one that finds the trailing space, because 'ok' and 'ok ' appear as two rows; add normalize=True and it gives proportions, which is how you notice that twenty percent of calls errored. And sample beats head, because the first five rows of a file are often a header artefact or the oldest and least representative data in it.

calls["status"].value_counts(normalize=True, dropna=False)
# ok       0.8
# error    0.2

calls.duplicated(subset=["request_id"]).sum()     # are the keys actually unique
calls["created_at"].agg(["min", "max"])           # is the date range what you think
calls.isna().mean().sort_values(ascending=False)  # fraction missing, worst first

Filter and derive#

ok = calls[calls["status"] == "ok"].copy()
ok["cost_usd"] = ok["input_tokens"] * 3e-6 + ok["output_tokens"] * 15e-6   # vectorised
ok["slow"] = ok["latency_ms"] > 5_000

Derive columns with arithmetic on whole columns, never with apply over rows. The difference is not marginal. Computing a product two ways on two hundred thousand rows: apply with axis=1 took 1,445 ms and the vectorised expression took 0.32 ms, a factor of about four and a half thousand, because apply builds a Python Series object for every row.

# 1,445 ms on 200,000 rows
df["c"] = df.apply(lambda r: r["a"] * r["b"], axis=1)

# 0.32 ms on the same rows
df["c"] = df["a"] * df["b"]

# For branching, np.where and np.select stay vectorised:
import numpy as np
ok["tier"] = np.where(ok["cost_usd"] > 0.01, "expensive", "cheap")
ok["band"] = np.select(
    [ok["latency_ms"] < 500, ok["latency_ms"] < 2000],
    ["fast", "normal"],
    default="slow",
)

Group and aggregate#

per_model = ok.groupby("model", observed=True).agg(
    calls=("user_id", "size"),
    users=("user_id", "nunique"),
    cost=("cost_usd", "sum"),
    p95_ms=("latency_ms", lambda s: s.quantile(0.95)),
).sort_values("cost", ascending=False)

#         calls  users     cost  p95_ms
# model
# sonnet      2      2  0.01290  817.00
# haiku       2      2  0.00225  154.25

Three details are worth knowing here. size counts rows and count counts non-null values in the named column, so on a column with nulls they disagree and the wrong one gives you a denominator that is too small. The group key becomes the index; add as_index=False, or call reset_index(), when you want it back as a column for a later merge. And to keep every original row while attaching a per-group figure, use transform rather than agg, which is how a share-of-group column is computed in one line.

ok["model_cost"] = ok.groupby("model", observed=True)["cost_usd"].transform("sum")
ok["share_of_model"] = ok["cost_usd"] / ok["model_cost"]

# Top 3 calls per user by cost, the pandas spelling of the SQL window pattern:
top3 = (ok.sort_values("cost_usd", ascending=False)
          .groupby("user_id", observed=True)
          .head(3))

Joins that check themselves#

users = pd.read_csv("users.csv")
merged = ok.merge(users, on="user_id", how="left", validate="many_to_one", indicator=True)
print(merged["_merge"].value_counts())      # how many calls found no user?
# both          4
# left_only     1
# right_only    0

validate raises an error if users has duplicate ids, which would otherwise silently multiply your rows and your cost totals. indicator shows how many rows failed to match. Without validate, five rows joined to a users table with one duplicate id come back as seven, every sum over them is inflated, and there is no error and no warning — only a revenue figure that is wrong in a way nobody notices for a month.

# MergeError: Merge keys are not unique in right dataset;
#             not a many-to-one merge

# When it fires, find them and decide what one row per key means:
dupes = users[users.duplicated("user_id", keep=False)].sort_values("user_id")
users_latest = (users.sort_values("updated_at")
                     .drop_duplicates("user_id", keep="last"))

Time#

daily = ok.set_index("created_at").resample("D")["cost_usd"].sum()
# 2026-03-01    0.00885
# 2026-03-02    0.00540
# 2026-03-03    0.00000     <- days with no calls are present, as zero
# ...
# 2026-03-09    0.00090

weekly_by_model = ok.pivot_table(index=pd.Grouper(key="created_at", freq="W"),
                                 columns="model", values="cost_usd", aggfunc="sum",
                                 observed=True)
# model         haiku  sonnet
# created_at
# 2026-03-01  0.00135  0.0075
# 2026-03-08      NaN  0.0054     <- NaN, not 0: no haiku calls that week

Two things to get right about time. resample needs a datetime index, which is what the set_index is for; Grouper does the same job inline when the frame has other grouping keys. And a frequency string means the period it ends on: W is weeks ending Sunday by default, ME is month end, MS is month start. Getting that wrong shifts every figure by a period, which is hard to see and easy to present.

Store timestamps in UTC and convert only when displaying. A daily total computed in local time counts the wrong hours on the days the clocks change, and comparing two days across a daylight saving boundary compares 23 hours with 25.

calls["created_at"] = pd.to_datetime(calls["created_at"], utc=True, errors="coerce")
calls["created_at"].isna().sum()     # rows whose date would not parse: check this

local = calls["created_at"].dt.tz_convert("Europe/London")
calls["hour"] = local.dt.hour
calls["date"] = local.dt.date

You can already do the work in pandas. This is what to do when a step that used to take seconds takes twenty minutes, or dies with a memory error. Four things usually fix it: replacing a Python loop with a vectorised expression, storing the data as Parquet rather than CSV, and moving to Polars or DuckDB, both of which run on one machine and handle files larger than memory. Knowing them means you reach for a cluster when you genuinely need one, which is far rarer than people assume.

Measure#

Find the slow step before changing anything, because the answer is routinely not where people assume. Time the stages separately rather than the whole script, and look at memory as well as seconds: a job that swaps to disk looks like a slow job and is really a memory problem.

%timeit df["total"] = [a * b for a, b in zip(df["qty"], df["price"])]   # Python loop
%timeit df["total"] = df["qty"] * df["price"]                            # vectorised

# Measured on 200,000 rows:
#   df.apply(..., axis=1)   1,445 ms
#   df["a"] * df["b"]           0.32 ms      about 4,500x

import time
t = time.perf_counter()
df = pd.read_parquet("calls.parquet")
print(f"load {time.perf_counter() - t:.1f}s  {df.memory_usage(deep=True).sum()/1e9:.1f} GB")

What is slow                          What actually fixes it
a per-row apply or an explicit loop   vectorise, or np.select for branching
reading the same CSV repeatedly       convert to Parquet once
reading 30 columns to use 3           Parquet plus a column list
merging on a string key               narrow the dtype, or merge on an integer id
it fits, but only just                narrower dtypes and categoricals
it does not fit at all                Polars lazy, or DuckDB over the files
none of the above                     then, and only then, a cluster

Parquet#

CSV is text: every read parses every value. Parquet stores typed columns, compressed, and lets readers load only the columns and row groups they need. Converting a large CSV to Parquet once often makes every later step several times faster and the files much smaller.

Three properties do the work. It is columnar, so reading two columns of thirty reads two columns of data rather than every byte of every row. It carries a schema, so a numeric column comes back numeric with no inference, and a type error in the data is caught at write time rather than discovered a week later. And it stores statistics per row group, so a reader can skip whole blocks that cannot match a filter. Partitioning the files by date on disk takes that further: a query for one day opens one directory.

df.to_parquet("calls.parquet", compression="zstd")
df.to_parquet("calls/", partition_cols=["date"])     # calls/date=2026-03-11/...

pd.read_parquet("calls.parquet", columns=["user_id", "cost_usd"])   # two columns only
pd.read_parquet("calls/", filters=[("date", "=", "2026-03-11")])    # one partition

Polars#

import polars as pl

result = (
    pl.scan_parquet("calls/*.parquet")            # lazy: nothing is read yet
      .filter(pl.col("status") == "ok")
      .group_by("model")
      .agg(pl.col("latency_ms").quantile(0.95).alias("p95"), pl.len().alias("calls"))
      .collect()                                  # plan optimised, then run in parallel
)

The word doing the work is scan rather than read. Nothing happens until collect, so Polars sees the whole chain first and can push the filter down into the file reader and load only the columns the query mentions. A pandas version of the same chain reads every column of every file into memory and then discards most of it, which is why the difference is often large rather than marginal. Polars also uses every core by default, where pandas uses one.

The API is different enough to be worth an hour of adjustment. There is no index, which removes a whole category of alignment surprises. Expressions are built with pl.col rather than by subscripting the frame, so a filter is written as a value passed to filter rather than a boolean frame. And there is no in-place mutation: you add columns with with_columns and get a new frame back.

pandas                                    Polars
df[df.status == "ok"]                     df.filter(pl.col("status") == "ok")
df["c"] = df["a"] * df["b"]               df.with_columns((pl.col("a") * pl.col("b")).alias("c"))
df.groupby("k").agg(n=("a", "sum"))       df.group_by("k").agg(pl.col("a").sum().alias("n"))
pd.read_parquet(p)                        pl.read_parquet(p), or scan_parquet for lazy
df.merge(o, on="k", how="left")           df.join(o, on="k", how="left")
df.to_pandas() is not needed              df.to_pandas() when a library demands pandas

DuckDB#

import duckdb

duckdb.sql("""
    SELECT model, count(*) AS calls, quantile_cont(latency_ms, 0.95) AS p95
    FROM 'calls/*.parquet'
    WHERE status = 'ok'
    GROUP BY model
""").df()

DuckDB is a database engine that runs inside your process with no server to install and no data to load. It reads Parquet, CSV and JSON files directly, and it also reads a pandas or Polars frame that is already in your session by name, which makes it easy to drop into an existing script for the one step that is slow. It streams, so it will happily group by over files larger than memory, and it uses every core.

calls = pd.read_parquet("calls.parquet")     # an ordinary pandas frame

duckdb.sql("SELECT model, sum(cost_usd) FROM calls GROUP BY model").df()
# 'calls' is found in the enclosing Python scope; nothing was copied in

duckdb.sql("COPY (SELECT * FROM 'raw/*.csv') TO 'out.parquet' (FORMAT parquet)")
# convert a directory of CSVs to Parquet without loading any of it

con = duckdb.connect("analysis.duckdb")      # or persist it to a file

Choosing between them is mostly about which language you want to write. DuckDB when the operation is naturally a query and especially when it joins several files, Polars when it is a chain of transformations inside Python code. Both run on one machine, and a single machine with fast disks handles tens of gigabytes comfortably, which is more than most jobs that get moved to a cluster.