Axes and Reductions

In one line: axis=k means "collapse axis k" — and getting axis=0 versus axis=1 backwards is the most consequential silent error in matrix bioinformatics.

The rule that actually works

Forget "rows" and "columns" for a moment. The axis you name is the axis that disappears.

a = np.arange(12).reshape(3, 4)     # shape (3, 4)

a.sum(axis=0).shape      # (4,)   — axis 0 (length 3) is gone
a.sum(axis=1).shape      # (3,)   — axis 1 (length 4) is gone
a.sum().shape            # ()     — everything gone, scalar
a.sum(axis=(0, 1))       # same as a.sum()

Translating back: for a 2-D array, axis=0 collapses down the rows and gives you one value per column; axis=1 collapses across the columns and gives one value per row.

Say it in terms of the output, not the input. "I want one value per gene, and genes are rows, so the output length is n_genes, so I collapse axis 1."

keepdims

a.mean(axis=1).shape                    # (3,)
a.mean(axis=1, keepdims=True).shape     # (3, 1)
a - a.mean(axis=1, keepdims=True)       # broadcasts cleanly

Use keepdims=True whenever the result feeds back into an operation with the original array. It removes the need for [:, None] and eliminates a whole class of Broadcasting errors.

The reduction catalogue

a.sum(); a.prod(); a.mean(); a.std(); a.var()
a.min(); a.max(); a.ptp()               # ptp = peak-to-peak (max - min)
a.argmin(); a.argmax()                  # index of the extreme
a.any(); a.all()
np.median(a, axis=1)
np.percentile(a, [25, 50, 75], axis=1)
np.quantile(a, 0.95, axis=1)
np.count_nonzero(a, axis=1)

a.cumsum(axis=1); a.cumprod(axis=1)     # not reductions — shape preserved
np.diff(a, axis=1)                       # shape shrinks by 1 along that axis

NaN-aware versions skip missing values instead of propagating them:

np.nansum, np.nanmean, np.nanstd, np.nanmedian,
np.nanmin, np.nanmax, np.nanpercentile, np.nanargmax

See Missing Data in NumPy.

The ddof trap

a.std()              # ddof=0 → population SD (divides by n)
a.std(ddof=1)        # sample SD (divides by n-1)

NumPy defaults to population (ddof=0). pandas defaults to sample (ddof=1). So np.std(x) and pd.Series(x).std() give different answers, and neither library is wrong. In statistics you almost always want ddof=1. Be explicit in both.

argmin / argmax and flat indices

a.argmax()                       # index into the FLATTENED array
np.unravel_index(a.argmax(), a.shape)    # → (row, col)
a.argmax(axis=1)                 # per-row column index

Ties go to the first occurrence. argmax on a boolean array finds the first True — but returns 0 when there is no True at all, so always guard with mask.any().

Bioinformatics: the two directions

With expr of shape (n_genes, n_samples):

expr.mean(axis=1)        # per-gene mean expression      → (n_genes,)
expr.sum(axis=0)         # per-sample library size       → (n_samples,)
expr.std(axis=1, ddof=1) # per-gene variability          → (n_genes,)
(expr > 0).sum(axis=1)   # in how many samples is each gene detected?
(expr > 0).sum(axis=0)   # how many genes detected per sample?
expr.argmax(axis=1)      # which sample expresses each gene most?

The sanity check: print the shape of every reduction result and confirm it equals the number of things you expect. (20000,) = genes. (500,) = samples. If you get 500 when you wanted 20000, you have the axis backwards.

This matters because with a square-ish matrix, or with a downstream step that broadcasts either way, the wrong axis produces a plausible result. Per-sample means and per-gene means are both smooth, reasonable-looking numbers. There is no error message.

Reductions along more than two axes

img = np.zeros((10, 256, 256))       # 10 z-slices of a microscopy image
img.mean(axis=0)                     # z-projection → (256, 256)
img.mean(axis=(1, 2))                # per-slice mean intensity → (10,)
img.max(axis=0)                      # maximum-intensity projection

Segmented and grouped reductions

When groups are not whole axes:

np.add.reduceat(values, segment_starts)      # sum within contiguous segments
np.bincount(group_codes, weights=values)     # sum by integer group id
np.bincount(group_codes)                      # count by group

For labelled groups, this is where you switch to pandas groupby. See Split-Apply-Combine.

Common mistakes

• axis=0 vs axis=1 reversed. Silent, plausible, wrong. Always check the output shape.
• Forgetting keepdims=True and then fighting broadcasting.
• ddof mismatch between NumPy and pandas in the same analysis.
• mean() on data containing NaN — the result is NaN, and one bad cell poisons the whole row. Use nanmean, but know that it silently changes your denominator.
• argmax on an all-False mask returning 0, which looks like a valid index.
• sum() on int32 counts overflowing. Pass dtype=np.int64.
• Assuming axis=-1 is "columns". It is the last axis, which is columns only for 2-D. It is the right choice when you want "along the innermost dimension" regardless of rank.

See also

Broadcasting · ufuncs · Missing Data in NumPy · ndarray · Split-Apply-Combine · groupby · Applied - Expression Matrices with NumPy