NumPy Mathematical Functions — The NaN Trap in np.mean

If you've ever run np.mean on a dataset with missing values and gotten nan back, you've hit the NaN trap. This default behavior silently propagates missing data through your pipeline, corrupting downstream aggregates and model predictions. The fix is simple: use np.nanmean or add defensive checks before aggregation. Ignoring this distinction will cost you hours of debugging production systems.

Why np.mean Is Not a Simple Average

NumPy mathematical functions are compiled vectorized operations that apply element-wise or aggregate computations over ndarrays without Python-level loops. The core mechanic is that these functions, including np.mean, operate on contiguous memory blocks via C/Fortran routines, achieving near-constant overhead per array regardless of size. This means np.mean does not compute a naive sum divided by count — it uses a numerically stable pairwise summation algorithm that reduces floating-point error.

In practice, np.mean accepts an axis argument for multi-dimensional reduction, a dtype parameter to control accumulator precision, and — critically — treats NaN values as propagators by default. If any element in the array is NaN, the result is NaN. This is not a bug; it's IEEE 754 floating-point semantics. The function does not silently skip missing data unless you explicitly call np.nanmean.

Use np.mean when you need a fast, numerically robust average over clean data. In real systems — sensor telemetry, financial tick data, ML feature pipelines — NaN values are common. Relying on np.mean without NaN handling leads to silent NaN propagation, corrupting downstream aggregates and model predictions. Always validate or sanitize inputs, or switch to np.nanmean for production pipelines.

Universal Functions (ufuncs)

Ufuncs apply a function to every element of an array. They support broadcasting, type casting, and output array specification.

import numpy as np

a = np.array([1.0, 4.0, 9.0, 16.0])

print(np.sqrt(a))          # [1. 2. 3. 4.]
print(np.log(a))           # [0.    1.386 2.197 2.773]
print(np.exp([0, 1, 2]))   # [1.    2.718 7.389]
print(np.abs([-3, -1, 2])) # [3 1 2]

# Two-array ufuncs
b = np.array([2.0, 3.0, 4.0, 5.0])
print(np.add(a, b))       # same as a + b
print(np.maximum(a, b))   # element-wise max
print(np.power(b, 2))     # [4. 9. 16. 25.]

Aggregation Functions and the axis Parameter

axis=0 operates down rows (along columns). axis=1 operates across columns (along rows). axis=None reduces everything to a scalar.

import numpy as np

m = np.array([[1, 2, 3],
              [4, 5, 6]])

print(m.sum())           # 21 — all elements
print(m.sum(axis=0))     # [5 7 9] — column sums
print(m.sum(axis=1))     # [6 15] — row sums

print(m.mean(axis=0))    # [2.5 3.5 4.5]
print(m.max(axis=1))     # [3 6]
print(m.argmax(axis=1))  # [2 2] — index of max in each row

Statistical Functions

NumPy covers standard statistics. For variance and standard deviation, note the ddof parameter — ddof=0 is population (default), ddof=1 is sample.

import numpy as np

data = np.array([2, 4, 4, 4, 5, 5, 7, 9])

print(np.mean(data))          # 5.0
print(np.median(data))        # 4.5
print(np.std(data))           # 2.0  (population)
print(np.std(data, ddof=1))   # 2.138 (sample)
print(np.var(data))           # 4.0
print(np.percentile(data, 75)) # 5.75
print(np.corrcoef([1,2,3], [1,2,3]))  # correlation matrix

Broadcasting: How ufuncs Handle Different Shapes

Broadcasting allows ufuncs to operate on arrays of different shapes by automatically expanding dimensions. The arrays must be compatible: they align from the rightmost dimension, and each dimension must be equal or one of them must be 1.

import numpy as np

a = np.array([[1, 2, 3],
              [4, 5, 6]])  # shape (2,3)
b = np.array([10, 20, 30]) # shape (3,) -> broadcasts to (1,3) then (2,3)
print(a + b)  # element-wise addition
# [[11 22 33]
#  [14 25 36]]

# Scalar broadcasting works too
print(a * 2)  # multiplies every element by 2

Performance: Vectorisation vs Python Loops

Vectorised operations in NumPy run in compiled C and are orders of magnitude faster than Python loops. The performance gap grows with array size — for a 10-million-element array, vectorised sum takes ~5ms, while a Python loop takes over a minute.

import numpy as np
import time

arr = np.random.randn(10_000_000)

# Vectorised
start = time.time()
result = np.sum(arr)
print(f"Vectorised: {time.time() - start:.5f}s")

# Python loop
start = time.time()
total = 0.0
for x in arr:
    total += x
print(f"Loop: {time.time() - start:.5f}s")

Advanced ufunc Methods: reduce, accumulate, outer

Ufuncs provide methods beyond direct element-wise application. reduce applies the operation cumulatively along an axis, accumulate returns all intermediate results, and outer computes the outer product.

import numpy as np

a = np.array([1, 2, 3, 4])

# reduce: cumulative sum produces single value
print(np.add.reduce(a))  # 10 = 1+2+3+4

# accumulate: all intermediate sums
print(np.add.accumulate(a))  # [1 3 6 10]

# outer: outer product
print(np.multiply.outer([1,2,3], [10,20,30]))
# [[10 20 30]
#  [20 40 60]
#  [30 60 90]]

Handling NaN in Statistical Aggregations

By default, NumPy aggregation functions (mean, sum, std) return NaN if any element is NaN. Use nan-aware versions: np.nansum, np.nanmean, np.nanstd. They ignore NaN values and compute on the remaining elements.

import numpy as np

data = np.array([1.0, 2.0, np.nan, 4.0])

print(np.mean(data))   # nan
print(np.nanmean(data)) # 2.333...

print(np.std(data))     # nan
print(np.nanstd(data))  # ~1.247

# For percentage of missing values
print(np.isnan(data).sum())  # 1

Trigonometric Functions: Radians or Regret

NumPy's trig functions – sin, cos, tan – take radians, not degrees. This is the #1 rookie mistake that costs hours in debugging orbital mechanics or signal processing. You pass degrees, you get garbage. No error, no warning. Just silent corruption.

Use np.deg2rad() or the shorthand np.radians() before feeding angles to sin/cos/tan. Inverse functions (arcsin, arccos, arctan) return radians. Convert back with np.rad2deg() if you need human-readable angles.

The arctan2(y, x) variant is your friend for full quadrant resolution. It handles the sign of both arguments, giving you the correct angle from -π to π. Single-argument arctan only returns -π/2 to π/2. That ambiguity has broken more than a few path-finding algorithms.

hypot computes Euclidean norm directly – use it instead of sqrt(x2 + y2). It avoids overflow for large values. unwrap fixes phase jumps in continuous signals by adding/subtracting multiples of 2π.

// io.thecodeforge — python tutorial

import numpy as np

# Probably wrong: degrees where radians expected
angles_deg = np.array([0, 30, 45, 60, 90])
wrong = np.sin(angles_deg)  # sin(30 deg) != sin(30 rad)
print(wrong)
# [ 0.        -0.98803162  0.85090352 -0.30481062  0.89399666]

# Correct: convert to radians first
angles_rad = np.radians(angles_deg)
right = np.sin(angles_rad)
print(right)
# [0.         0.5        0.70710678 0.8660254  1.        ]

# Full quadrant with arctan2
points = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]])
correct_angles = np.arctan2(points[:, 1], points[:, 0])
print(np.degrees(correct_angles))
# [ 45. 135. -135. -45.]

# Hypotenuse without overflow
leg_a, leg_b = 1e200, 1e200
print(np.hypot(leg_a, leg_b))  # safe
# 1.4142135623730951e+200

Rounding Functions: Know Which Floor You're On

Rounding looks trivial until you need precision control. NumPy gives you several options, and picking the wrong one will silently truncate your data.

np.round() rounds to the nearest even number – banker's rounding. It's the default. Use it for financial calculations where tie-breaking bias matters. But don't assume it rounds "up" for .5; it rounds to even. 2.5 becomes 2, 3.5 becomes 4.

np.floor() always goes down. np.ceil() always goes up. These are deterministic and predictable. Use them when you need a ceiling for bin edges or a floor for price thresholds.

np.trunc() just chops off the decimal part – equivalent to int() but vectorised. It's the fastest option if you only need integer parts.

np.rint() rounds to the nearest integer but returns a float array. Useful when you need the number type to stay consistent. np.fix() does the same as trunc but handles negative values correctly – it rounds toward zero.

Don't confuse np.around with np.round – they're the same function. And don't use Python's built-in round() on arrays; it doesn't vectorise.

// io.thecodeforge — python tutorial

import numpy as np

values = np.array([2.5, 3.5, -2.5, -3.5, 2.51])

print("Original:", values)
print("round (bankers):", np.round(values))
print("floor:", np.floor(values))
print("ceil:", np.ceil(values))
print("trunc:", np.trunc(values))
print("rint:", np.rint(values))
print("fix:", np.fix(values))

# Production example: binning sensor data
raw_temps = np.array([23.7, 24.1, 25.0, 25.9])
# Bin to nearest 0.5 degree (round to one decimal, then scale)
binned = np.round(raw_temps * 2) / 2
print("Binned temps:", binned)

Exponentials and Logarithms: Watch the Domain

Exponential and log functions are cheap in NumPy until you hit an edge case. Then they silently return inf or -inf, and your pipeline produces garbage.

np.exp(x) computes e^x element-wise. For x > 709, you overflow to infinity. For x < -745, you underflow to zero. Neither crashes – both silently produce values that will ruin downstream calculations. Check your input range before calling exp.

np.log(x) computes natural log. Domain: x > 0. Pass zero or negative and you get -inf or nan. log of zero happens more often than you think in real-world data. np.log1p(x) computes log(1+x) accurately for small x. Use it when x could be near zero or negative.

np.log10 and np.log2 exist for base-specific needs. They're slightly faster than dividing by log(10) or log(2).

np.expm1(x) computes e^x - 1 accurately for small x. Pair it with log1p for financial calculations where tiny differences matter.

Real production trap: power data. If you compute np.exp(-5 * t) for a range of t, the values decay to zero fast. But np.log(normalized_data) on data that's been quantized to zero kills your variance.

// io.thecodeforge — python tutorial

import numpy as np

# Overflow: gone silent large
large = np.array([700, 710, 800])
print(np.exp(large))
# [ 1.01423205e+304  1.51932128e+308                inf]

# Underflow: gone silent small
small = np.array([-745, -800, -1000])
print(np.exp(small))
# [9.88131292e-324  0.00000000e+000  0.00000000e+000]

# Log domain errors
bads = np.array([-1, 0, 1e-10])
print(np.log(bads))
# [ nan -inf -23.02585093]

# log1p fixes small x precision
x = np.array([1e-10, 1e-20, 1e-30])
print("log(1+x):", np.log(1 + x))
print("log1p:", np.log1p(x))

# expm1 for small x
print("exp(x)-1:", np.exp(x) - 1)
print("expm1:", np.expm1(x))

Hyperbolic Sine: Not Your High School Trig

Forget SOH CAH TOA. Hyperbolic sine (sinh) models real-world physical systems: hanging cables, heat transfer, special relativity. Unlike circular sine which oscillates, sinh grows exponentially. That means small inputs stay tame, but blow up fast past x=10 — overflow territory.

NumPy's np.sinh runs as a vectorized ufunc. It expects radians (hyperbolic ones have no unit, but numpy doesn't care about your confusion). Production gotcha: tanh is bounded [-1,1]; sinh is unbounded. Use np.clip before plotting or feeding into neural networks to prevent NaN cascades from massive values. Watch your domain — sinh(0) is 0, sinh(1) is ~1.175, sinh(20) is ~2.4e8.

// io.thecodeforge — python tutorial

import numpy as np

x = np.array([0, 1, 2, 10, 20], dtype=np.float64)
y = np.sinh(x)

print("sinh:", y)
print("Overflow risk?", np.any(np.isinf(y)))

# Clamp before using in further computation
y_safe = np.clip(y, -1e6, 1e6)
print("Clipped:", y_safe)

Cube Root: The Misunderstood Power

Square roots get all the fame. Cube roots handle negative numbers correctly — no imaginary units, no domain errors. np.cbrt(x) is faster and more numerically stable than x(1/3). Why? (1/3) risks floating-point errors with negative bases (Python returns a complex number). cbrt stays real and vectorized.

Use it for scaling features with heavy skew, volumetric calculations, or normalizing data with cubic relationships. Production advice: if you manually compute x**(1/3), you're writing a bug for negative inputs. Switch to np.cbrt — it's a dedicated ufunc with better precision. For integer arrays, cbrt returns float64 by default; cast back if needed.

// io.thecodeforge — python tutorial

import numpy as np

vals = np.array([-27, -8, 0, 8, 27], dtype=np.float64)
cubes = np.cbrt(vals)

print("Cube roots:", cubes)

# Compare with power
bad = vals ** (1/3)
print("Power roots:", bad)

# Notice: 1/3 is not exact in float, small errors visible
print("Difference:", cubes - bad)

Linear Algebra: Thinking in Matrices, Not Loops

NumPy's linear algebra module (np.linalg) moves you from element-wise operations to whole-matrix mathematics. Why this matters: Most data science pipelines — regression, PCA, neural nets — are linear algebra under the hood. Writing matrix multiplication with loops is bug-prone and 50x slower. np.dot or @ operator computes dot products, matrix-vector products, and matrix-matrix products in optimized C. np.linalg.inv inverses matrices; np.linalg.eig finds eigenvalues. Common pitfall: multiplying two 2D arrays with * does element-wise multiplication, not matrix multiplication. Use @. Avoid explicit loops for any linear transform. np.linalg.solve solves linear systems Ax=b directly, replacing manual Gaussian elimination. Always check matrix shape compatibility before multiplication.

// io.thecodeforge — python tutorial

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Matrix multiplication (not element-wise)
matmul = A @ B  # or np.dot(A, B)
print(matmul)
# Solve linear system Ax = b
b = np.array([9, 10])
x = np.linalg.solve(A, b)
print(x)
# Eigenvalues
eigvals = np.linalg.eigvals(A)
print(eigvals)

Matrix Multiplication: Why @ Beats Nested Loops

Matrix multiplication is the backbone of neural networks and linear models. In pure Python, a 1000x1000 matrix multiply takes nested loops (O(n³)) and minutes. NumPy's @ operator or matmul() runs compiled BLAS libraries — 100x faster. Why this matters: you can't scale models with manual loops. The key insight: matrix multiplication requires the inner dimensions to match (e.g., (3,4) @ (4,2) yields (3,2)). Shape errors are silent until runtime — validate with .shape. Avoid np.dot for new code; @ is cleaner. For batched multiplication (e.g., 1000 images at once), use 3D arrays with @ — it broadcasts across the batch dimension automatically. Production tip: if matrices exceed RAM, use np.memmap with chunked multiplication.

// io.thecodeforge — python tutorial

import numpy as np

# Shape (2,3) @ (3,2) -> (2,2)
X = np.array([[1, 2, 3], [4, 5, 6]])
Y = np.array([[7, 8], [9, 10], [11, 12]])

result = X @ Y
print(result)
print(result.shape)

# Batched: (2, 2, 3) @ (2, 3, 2) -> (2, 2, 2)
batch_X = np.stack([X, X])
batch_Y = np.stack([Y, Y])
batch_result = batch_X @ batch_Y
print(batch_result.shape)

Basic Arithmetic Operations: Warp Speed Math

Forget Python's scalars — NumPy arrays perform element-wise arithmetic without loops, and that's the whole point. Adding two arrays with + or np.add triggers the same compiled C loops that make NumPy fast. Subtraction, multiplication, and division work identically, broadcasting smaller arrays into compatible shapes automatically. The why is performance: a 10-million-element addition using + runs in milliseconds; a Python for-loop takes seconds. np.divide handles integer division carefully, returning floats by default (unlike Python's //). Use np.floor_divide when you want the floor after division. For mod, np.mod or % gives the remainder, but watch for negative dividends — NumPy follows C's convention (sign of divisor). np.fmod follows C's fmod (sign of dividend). These operators are not syntactic sugar; they are the backbone of vectorised computation. Avoid Python's for lists (repetition) — NumPy arrays don't repeat, they multiply element-wise. Dot product? That's @ or np.dot, not .

// io.thecodeforge — python tutorial
// 25 lines max
import numpy as np

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

# Element-wise operations — no loops
sum_arr = a + b          # [5, 7, 9]
diff_arr = np.subtract(a, b)  # [-3, -3, -3]
prod_arr = a * b         # [4, 10, 18]
div_arr = a / b          # [0.25, 0.4, 0.5]

# Integer division: floor vs truncation
floor_div = np.floor_divide(a, b)  # [0, 0, 0]

# Modulo: careful with negatives
c = np.array([-1, -2, -3])
print(c % 2)   # [1, 0, 1] (sign of divisor)
print(np.fmod(c, 2))  # [-1, 0, -1] (sign of dividend)

Exponentiation: Power Up, Not Loop Down

Raising numbers to powers in plain Python uses , but NumPy's np.power and operator are vectorised and accelerated. np.power(base, exp) computes base^exp element-wise, with broadcasting support for scalar exponents. Why use np.power over ? For integer exponents, np.power allows dtype control: np.power(x, 2, dtype=np.float64) avoids overflow on large ints. For fractional exponents (e.g., square root), use np.sqrt (faster than 0.5). np.exp computes e^x, useful for softmax and logistic regression — but that's exponential (section not covered here). Exponentiation with arrays as exponents is common in physics: np.power(x, y) where both are arrays. Beware of 0^0 returns 1 (consistent with IEEE 754). Negative base with fractional exponent? NumPy returns nan for real arrays; use complex dtype. For speed, np.exp uses the exp C library call, not a Python loop. Avoid writing x**2 for large arrays where you can use np.square(x) — it's slightly faster by skipping general power logic.

// io.thecodeforge — python tutorial
// 25 lines max
import numpy as np

x = np.array([1, 2, 3, 4])

# Power with scalar
squares = x ** 2                # [1, 4, 9, 16]
cubes = np.power(x, 3)          # [1, 8, 27, 64]

# Power with array exponent
exps = np.array([0, 1, 2, 3])
result = np.power(x, exps)      # [1^0, 2^1, 3^2, 4^3] = [1, 2, 9, 64]

# Use np.square for speed
fast_sq = np.square(x)          # same as x**2

# Fractional: use sqrt or power
sqrt_via_power = np.power(x, 0.5)
sqrt_via_sqrt = np.sqrt(x)      # faster
print(sqrt_via_sqrt)