Linear Algebra with NumPy

Core Principle: High-performance linear algebra operations in NumPy execute via the @ operator and np.linalg routines, backed by optimized C and Fortran BLAS/LAPACK implementations.

Multiplication: three different things

a * b                # ELEMENTWISE: this is not matrix multiplication
a @ b                # matrix multiplication
np.matmul(a, b)      # same as @
np.dot(a, b)         # same for 2-D; different for >2-D: avoid
a.dot(b)             # method form
np.einsum("ij,jk->ik", a, b)   # explicit index notation

* is elementwise. This is the number one confusion for people arriving from MATLAB or R, where * on matrices means matrix product. In NumPy, * broadcasts elementwise and @ does the linear algebra.

@ was added in Python 3.5 specifically for this. Use it: np.dot has subtle and surprising behaviour on 3-D+ inputs.

The core operations

from numpy import linalg as la

la.inv(a)                  # inverse: rarely what you actually want
la.pinv(a)                 # pseudo-inverse, handles singular matrices
la.solve(A, b)             # solve Ax = b: use this instead of inv(A) @ b
la.lstsq(A, b, rcond=None) # least squares for overdetermined systems
la.det(a)
la.matrix_rank(a)
la.norm(a); la.norm(a, axis=1)     # vector/matrix norms
la.eig(a); la.eigh(a)              # eigen (eigh for symmetric: faster, real)
la.svd(a, full_matrices=False)     # singular value decomposition
la.cholesky(a)                     # for positive-definite matrices
la.qr(a)

Never invert to solve. la.solve(A, b) is faster and numerically far more stable than la.inv(A) @ b. The inverse is almost never needed explicitly: if you think you need it, you probably need solve, lstsq, or pinv.

eigh over eig for symmetric/Hermitian matrices (covariance and correlation matrices are symmetric). It is faster and returns real, sorted eigenvalues instead of complex ones with numerical noise.

Correlation and covariance

np.corrcoef(m)                # rows are variables by default
np.corrcoef(m, rowvar=False)  # columns are variables
np.cov(m, rowvar=False, ddof=1)

The rowvar default trips people up constantly. For an expression matrix of shape (n_genes, n_samples):

np.corrcoef(expr)                      # gene–gene correlation → (n_genes, n_genes)
np.corrcoef(expr, rowvar=False)        # sample–sample        → (n_samples, n_samples)
np.corrcoef(expr.T)                    # same thing

A 20,000-gene correlation matrix is 20,000² × 8 bytes = 3.2 GB. Check your shapes before running this. Sample–sample correlation (500 × 500) is the one you usually want for QC.

PCA from scratch

The canonical dimensionality-reduction workflow, and a good test of whether you have internalised Axes and Reductions and Broadcasting:

# X: (n_samples, n_features): note the orientation!
X = log_expr.T                                  # samples as rows
Xc = X - X.mean(axis=0, keepdims=True)          # centre each feature

U, S, Vt = np.linalg.svd(Xc, full_matrices=False)

scores = U * S                    # (n_samples, n_components): plot these
loadings = Vt                     # (n_components, n_features): gene weights
var_explained = S**2 / np.sum(S**2)

print(f"PC1: {var_explained[0]:.1%}, PC2: {var_explained[1]:.1%}")

Two things people get wrong:

• Orientation. sklearn and most PCA conventions want samples as rows. Bulk expression matrices are stored genes-as-rows. Transpose first.
• Centring. SVD without centring gives you a first component that is just the overall mean. Always subtract the column means. Whether to also scale by SD (i.e. use correlation rather than covariance) depends on whether your features are on comparable scales: for log-CPM they usually are, so centring alone is standard.

Report the variance explained on your axis labels. A PCA plot without it is uninterpretable.

einsum

For operations awkward to express otherwise:

np.einsum("ij,jk->ik", a, b)      # matrix multiply
np.einsum("ii->i", a)              # diagonal
np.einsum("ij->j", a)              # column sums
np.einsum("ij,ij->i", a, b)        # row-wise dot products
np.einsum("ij,j->i", a, v)         # matrix-vector

Read the subscripts as: name each axis, repeat a name to sum over it, and list what survives after ->. The last example: row-wise dot products: is genuinely awkward without einsum ((a * b).sum(axis=1) works but materialises a temporary).

np.einsum(..., optimize=True) reorders multi-operand contractions to minimise work.

Performance notes

np.show_config() tells you which BLAS you are linked against (OpenBLAS, MKL, Accelerate). Matrix multiplication is multi-threaded there, which means:

• A single @ on large matrices already uses all your cores.
• Wrapping it in multiprocessing will oversubscribe and get slower. Set OMP_NUM_THREADS=1 in workers.
• float32 is roughly 2× faster than float64 for large matmuls and usually precise enough for genomics.

Bioinformatics uses

# sample-sample distance for a QC heatmap
from scipy.spatial.distance import pdist, squareform
d = squareform(pdist(log_expr.T, metric="correlation"))

# projecting new samples onto existing PCs
new_scores = (new_X - train_mean) @ Vt.T

# position weight matrix scoring: (4, L) PWM against one-hot windows (n, L, 4)
scores = np.einsum("nlb,bl->n", windows_onehot, pwm)

# linear model per gene, all genes at once (design matrix D, (n_samples, p))
betas = np.linalg.lstsq(D, log_expr.T, rcond=None)[0]     # (p, n_genes)

That last one is worth noting: lstsq solves for many right-hand sides simultaneously, so you can fit 20,000 per-gene linear models in one call. That is the computational core of limma.

Common mistakes

• * when you meant @. With square matrices this runs without error and gives the wrong answer.
• inv(A) @ b instead of solve(A, b).
• np.corrcoef orientation: 3.2 GB surprise.
• Not centring before PCA.
• Wrong PCA orientation (genes as observations).
• eig on a symmetric matrix, getting tiny imaginary components and confusion. Use eigh.
• Oversubscribing threads by parallelising over an already-threaded BLAS.
• Ignoring conditioning. la.cond(A) above ~1e10 means your solution is numerically unreliable regardless of the method.

See also

Axes and Reductions · Broadcasting · ndarray · Memory Layout and Strides · Applied - Expression Matrices with NumPy · Applied - Heatmaps and Clustermaps

Exercise

Write solve_system(a, b), which solves a @ x = b by whichever routine the matrix actually supports. Use np.linalg.solve only when a is square and of full rank; a rectangular or rank-deficient a must go to np.linalg.lstsq with rcond=None, because solve raises on a singular matrix rather than returning a useful answer.

Return four lines: method=<solve|lstsq> rank=<n>, using np.linalg.matrix_rank for the rank; x=<list> with each component rounded to three decimals; residual=<value>, the Euclidean norm of a @ x - b rounded to six decimals, which is zero for an exact solution and positive for a least-squares fit; and det=<value> rounded to three decimals, or det=NA when a is not square and no determinant exists.