Linear Algebra with NumPy

Analogy categoryAnalogies

NumPy arrays double as the natural representation of vectors and matrices, and the numpy.linalg submodule wraps battle-tested LAPACK/BLAS routines to perform the operations that make linear algebra practical at scale: matrix multiplication, inversion, determinants, solving linear systems, and eigen-decomposition. Understanding these operations matters because they underpin almost every quantitative technique used downstream in data science, from ordinary least squares regression to PCA and the forward pass of a neural network. Crucially, NumPy distinguishes between element-wise multiplication (the * operator) and true matrix multiplication (@ or np.matmul), a distinction that trips up many newcomers and silently produces wrong results if ignored.

🏏

Cricket analogy: NumPy's linalg submodule is like a stats analyst's toolkit for a season's run-and-wicket matrices, and confusing element-wise * with true matrix multiplication @ is like multiplying two batters' scores position-by-position instead of properly combining a lineup matrix with a bowling-attack matrix.

Matrix Multiplication and the @ Operator

Since Python 3.5, the @ operator is the idiomatic way to perform matrix multiplication on NumPy arrays, following the standard rule that an (m x n) matrix multiplied by an (n x p) matrix yields an (m x p) matrix, with the inner dimensions n required to match. np.matmul and np.dot behave the same way for 2-D arrays, but np.dot falls back to a scalar dot product for 1-D arrays, and np.matmul has slightly different broadcasting semantics for arrays with more than two dimensions (it treats extra leading dimensions as stacks of matrices). Using * on two 2-D arrays instead performs element-wise multiplication and requires the shapes to be broadcast-compatible rather than matrix-compatible, so mixing the two up is a common source of shape-mismatch bugs or, worse, silently wrong numeric results when shapes happen to broadcast.

🏏

Cricket analogy: The @ operator combines an (overs x players) scoring matrix with a (players x matches) weighting matrix into an (overs x matches) result, requiring the inner player counts to match, while plain * would just multiply corresponding cells, a common source of wrong totals if you mix them up.

python

import numpy as np

A = np.array([[1, 2], [3, 4], [5, 6]])   # shape (3, 2)
B = np.array([[2, 0], [1, 3]])           # shape (2, 2)

product = A @ B                          # matrix multiply -> shape (3, 2)
print(product)
# [[ 4  6]
#  [10 12]
#  [16 18]]

v = np.array([1, 0, -1])
print(np.dot(v, v))   # scalar dot product -> 2
print(A.T)            # transpose -> shape (2, 3)

Solving Linear Systems, Determinants, and Inverses

Given a system Ax = b, np.linalg.solve(A, b) computes x directly using an LU decomposition, and is both faster and more numerically stable than explicitly forming the inverse with np.linalg.inv(A) and multiplying. np.linalg.det(A) returns the determinant, which is zero (or numerically near-zero) exactly when A is singular and has no unique inverse; np.linalg.inv will raise a LinAlgError on an exactly singular matrix but can return a wildly inaccurate result on a near-singular one, so checking the determinant or condition number (np.linalg.cond) before trusting an inverse is good practice. For non-square or ill-conditioned systems, np.linalg.lstsq returns the least-squares solution instead.

🏏

Cricket analogy: np.linalg.solve() directly finds a batter's required run rate from a system of over-by-over equations, faster and more stable than explicitly inverting the matrix, while a near-zero determinant would flag an unsolvable, contradictory scorecard.

python

import numpy as np

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])

x = np.linalg.solve(A, b)
print(x)          # [2. 3.]  (solves 3x+y=9, x+2y=8)

print(np.linalg.det(A))   # 5.0 -> nonsingular, unique solution exists
print(np.linalg.inv(A))   # explicit inverse, generally avoid for solving systems

Eigenvalues, Eigenvectors, and Decompositions

np.linalg.eig(A) returns a tuple of eigenvalues and the matrix of corresponding eigenvectors (as columns), which is the foundation of techniques like Principal Component Analysis, where the eigenvectors of a covariance matrix identify directions of maximal variance. For symmetric or Hermitian matrices, np.linalg.eigh should be preferred: it exploits the guaranteed real eigenvalues and orthogonal eigenvectors of symmetric matrices to be both faster and numerically more stable than the general-purpose eig. NumPy also exposes np.linalg.svd for singular value decomposition, which generalizes eigen-decomposition to non-square matrices and underlies dimensionality reduction, pseudo-inverse computation, and recommender-system factorization.

🏏

Cricket analogy: np.linalg.eig() extracts the eigenvectors of a covariance matrix of batting stats, identifying the direction of maximal variance in a batter's performance, and np.linalg.eigh should be used since the covariance matrix is symmetric, giving faster, more stable results.

A useful mental model: @ is to matrices what * is to scalars — it composes linear transformations. Applying A @ B to a vector v first applies B, then A, exactly mirroring function composition f(g(v)).

A very common bug: writing A * B when you meant matrix multiplication. For two square matrices of the same shape, * will run without error (broadcasting kicks in silently) but produce element-wise products, not a matrix product — the code looks correct and even returns the 'right shape', but the numbers are wrong.

Use @ or np.matmul for true matrix multiplication; * is always element-wise.
np.linalg.solve(A, b) is preferred over computing inv(A) @ b for solving linear systems — faster and more numerically stable.
A zero (or near-zero) determinant signals a singular or ill-conditioned matrix; check np.linalg.cond before trusting an inverse.
Use np.linalg.eigh instead of np.linalg.eig whenever the matrix is symmetric, for speed and numerical stability.
Singular value decomposition (np.linalg.svd) generalizes eigen-decomposition to non-square matrices and underlies PCA and pseudo-inverses.
1-D arrays passed to np.dot compute a scalar dot product, not a matrix product.

Practice what you learned

1. What does the expression `A * B` compute for two 2-D NumPy arrays of the same shape?

2. Why is `np.linalg.solve(A, b)` generally preferred over `np.linalg.inv(A) @ b`?

3. What does a determinant of exactly zero indicate about a square matrix?

4. Which function should be preferred for eigen-decomposition of a symmetric matrix?

5. For 1-D arrays u and v, what does `np.dot(u, v)` return?

Was this page helpful?

Aggregations and Statistics

Topics covered

#Python#PandasNumPyStudyNotes#DataScience#LinearAlgebraWithNumPy#Linear#Algebra#NumPy#Matrix#StudyNotes#SkillVeris

Linear Algebra with NumPy