Axes and Reductions

Core Principle: Reductions aggregate values along a designated dimension (collapsing axis $k$), where axis=0 reduces across rows to summarize columns and axis=1 reduces across columns to summarize rows.

The rule that actually works

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

np.int64(66)

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

array([[-1.5, -0.5,  0.5,  1.5],
       [-1.5, -0.5,  0.5,  1.5],
       [-1.5, -0.5,  0.5,  1.5]])

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:

(<function nanmin at 0x1063a43b0>, <function nanmax at 0x1063a4530>, <function nanpercentile at 0x1063a4f30>, <function nanargmax at 0x1063a47f0>)

See Missing Data in NumPy.

The ddof trap

np.float64(3.605551275463989)

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

array([3, 3, 3])

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

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], shape=(256, 256))

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

Exercise

Write group_profile(matrix, labels). The matrix holds one gene per row and one sample per column; labels gives the group each column belongs to. Return ["genes=0"] when the matrix has no rows, which is the case a reduction over an empty axis would otherwise turn into a warning and a NaN.

Otherwise the first line is groups=<names joined by commas> in sorted order. Then emit one line per gene, row<i> <group> <mean>, naming the group whose columns give that gene its highest mean and that mean rounded to three decimals; where two groups tie, the earlier name in sorted order wins. Follow those with peak=row<i>,<group>, locating the single largest entry of the gene-by-group mean table by converting the flat argmax into a pair of indices, and finally spread=<max minus min> over that same table, rounded to three decimals.