NumPy: Math, Linear Algebra & Performance

Universal Functions (ufuncs)

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

# Element-wise math (vectorized — no Python loop needed)
np.sqrt(a)              # [1. 2. 3. 4.]
np.square(a)            # [1. 16. 81. 256.]
np.abs(np.array([-1, -2, 3]))  # [1 2 3]
np.exp(a)               # e^x
np.log(a)               # natural log
np.log2(a)
np.log10(a)
np.sin(a)
np.cos(a)
np.floor(a)
np.ceil(a)
np.round(a, decimals=2)
np.clip(a, 0, 5)        # clamp values to [0, 5]

# Aggregations
a.sum()                  # sum all elements
a.sum(axis=0)            # sum along rows (result: shape per column)
a.sum(axis=1)            # sum along columns (result: shape per row)
a.mean()
a.std()
a.var()
a.min()
a.max()
a.argmin()               # index of minimum
a.argmax()               # index of maximum
a.cumsum()               # cumulative sum
np.median(a)
np.percentile(a, [25, 50, 75])
np.unique(a)
np.sort(a)               # returns sorted copy
a.sort()                 # sorts in-place

Linear Algebra

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

# Matrix multiplication
A @ B                    # preferred (PEP 465)
np.matmul(A, B)          # same
np.dot(A, B)             # also works for 2D

# Element-wise (NOT matrix multiply)
A * B                    # [[ 5,12],[21,32]]

# Linear algebra functions
np.linalg.det(A)         # determinant
np.linalg.inv(A)         # inverse
np.linalg.trace(A)       # sum of diagonal
np.linalg.norm(A)        # Frobenius norm by default
np.linalg.norm(A, axis=1)  # row norms

# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)

# Singular Value Decomposition
U, S, Vt = np.linalg.svd(A)

# Solve linear system Ax = b
b = np.array([1, 2])
x = np.linalg.solve(A, b)     # faster than inv(A) @ b

# Least squares (overdetermined system)
A_tall = np.random.rand(100, 3)
b_tall = np.random.rand(100)
x, residuals, rank, sv = np.linalg.lstsq(A_tall, b_tall, rcond=None)

Performance & Best Practices

• Vectorize: never write Python loops over array elements. Use ufuncs, array operations, and built-in aggregations.

• Avoid copies: slices return views; fancy indexing and boolean indexing return copies. Profile with a.base to check.

• Memory layout: C-contiguous (row-major) for row operations; F-contiguous for column operations. Use np.ascontiguousarray() if needed.

• dtype matters: float32 is 2× faster than float64 on GPU and many CPU SIMD operations. Use smallest dtype that fits.

• np.einsum: Einstein summation notation for complex tensor operations — often faster than matmul chains.

• numba: JIT-compiles Python/NumPy to LLVM — near-C speed for loops that can't be vectorized.

• Use out= parameter to write results into pre-allocated array: np.add(a, b, out=result) — avoids allocation.

# Benchmark: Python loop vs NumPy vectorization
import time
n = 1_000_000
a = np.random.rand(n)

# Python loop — slow
start = time.time()
result = [x ** 2 for x in a]
print(f"Loop: {time.time() - start:.3f}s")  # ~0.3s

# NumPy vectorized — fast
start = time.time()
result = a ** 2
print(f"NumPy: {time.time() - start:.4f}s") # ~0.002s — 150× faster

# np.einsum examples
A = np.random.rand(100, 50)
B = np.random.rand(50, 30)
np.einsum('ij,jk->ik', A, B)   # matrix multiply (same as A @ B)
np.einsum('ij->i', A)          # row sums (same as A.sum(axis=1))
np.einsum('ii->', A[:50, :50]) # trace

# Save/load arrays
np.save('array.npy', a)                   # single array
np.savez('arrays.npz', a=a, b=B)         # multiple arrays
loaded = np.load('arrays.npz')
loaded['a']                               # retrieve by key
np.savetxt('data.csv', A, delimiter=',', fmt='%.4f')