Matrix multiplication, solving systems, and why you should almost never invert a matrix.
Overview
@ versus *
Linear algebra operations in NumPy trip people up because the notation looks almost right but means something different. The asterisk on two matrices gives elementwise multiplication. The at-sign gives matrix multiplication. Both produce valid arrays of plausible shape. Only one is the operation you intended.
Matrix multiplication has its own operator because it is fundamentally different from elementwise operations. Elementwise multiplication pairs corresponding elements and gives a result of the same shape. Matrix multiplication sums products along the inner dimension and gives a result whose shape depends on both inputs. For two (3, 3) arrays, both operations return (3, 3), but the values are entirely different.
The at-sign operator was added in Python 3.5 specifically for this purpose. Before that, you had to write np.dot or np.matmul, which are less readable in expressions with multiple terms. The three forms do the same thing for 2-D arrays. For higher dimensions, matmul and @ broadcast over batch dimensions, while dot does not.
Beyond the product, NumPy provides the common linear algebra operations in the np.linalg module. Solving a linear system is np.linalg.solve, which is faster and more numerically stable than computing the inverse and multiplying. Computing the inverse is np.linalg.inv. Determinant is np.linalg.det. Eigenvalues and eigenvectors are np.linalg.eig or np.linalg.eigh for symmetric matrices.
Checking the result of a linear algebra computation is important because floating-point arithmetic introduces rounding errors. The product of a matrix and its inverse is not exactly the identity matrix. It is close to it. np.allclose checks whether two arrays are equal within a tolerance, and it is the right way to verify such results.
This module covers the @ operator, np.dot, np.linalg.solve, inv, det, eig, and the norm function. The examples use small matrices where you can verify the result by hand. The focus is on the operations that appear in everyday data work, not on the full scope of numerical linear algebra.
* is elementwise. @ is mat product.
Both are valid on two square mats of the same size. They give fully other answers. Nothing raises, nothing warns. The result of the wrong one is a plausible-looking mat of the right shape.
This is the most common linear algebra bug in NumPy code. It lives review easily because both lines look correct.
np.matmul is the function form of @. np.dot agrees with it for 2-D. They diverge for higher dims. matmul broadcasts the leading axes and treats the last two as mats, dot does a thing other and rarely what you want. Prefer @ or matmul.
Worth knowing
@ is matrix multiplication; * is elementwise. Both are legal on the same arrays, so nothing warns you when you pick the wrong one.
Shapes must agree in the middle: (m,k) @ (k,n) gives (m,n).
A 1-D array is promoted to whichever vector orientation makes the multiplication valid.
To solve Ax = b, use np.linalg.solve — never inv(A) @ b. It is faster and numerically better.
matrix_rank and cond tell you whether a system is solvable; testing det == 0 in floating point does not.
norm for lengths and distances, eigvalsh for symmetric matrices, lstsq for fitting.
Linear Algebra
mat product, solving systems, and why you should almost never invert a mat.
@ is matrix multiplication, * is not
The single most common confusion, and both are legal so nothing warns you.
example_01.pyNumPy
import numpy as np
a = np.array([[1, 2],
[3, 4]])
b = np.array([[5, 6],
[7, 8]])
print("a * b (elementwise):"); print(a * b)
print()
print("a @ b (matrix product):"); print(a @ b)
print()
print("np.matmul and np.dot agree with @ for 2-D:")
print(" matmul same:", np.array_equal(a @ b, np.matmul(a, b)))
print()
print("Order matters - matrix multiplication is not commutative:")
print(" a@b == b@a :", np.array_equal(a @ b, b @ a))
Output
Shapes must line up in the middle
(m,k) @ (k,n) -> (m,n). The inner dimensions cancel.
example_02.pyNumPy
import numpy as np
a = np.ones((2, 3))
b = np.ones((3, 4))
print("(2,3) @ (3,4) ->", (a @ b).shape, " the 3s cancel")
try:
a @ np.ones((2, 4))
except ValueError as e:
print("(2,3) @ (2,4) -> ValueError")
print(" ", str(e)[:88])
print()
print("A 1-D array is treated as whichever vector makes the shapes work:")
v = np.array([1, 1, 1])
print(" (2,3) @ (3,) ->", (a @ v).shape, "matrix times column vector")
print(" (3,) @ (3,4) ->", (v @ b).shape, "row vector times matrix")
Output
Solving a system: A x = b
solve is the answer. Not inv.
example_03.pyNumPy
import numpy as np
# 2x + y = 11
# x + 3y = 18
A = np.array([[2.0, 1.0],
[1.0, 3.0]])
b = np.array([11.0, 18.0])
x = np.linalg.solve(A, b)
print("solution x, y :", np.round(x, 6))
print("check A @ x :", np.round(A @ x, 6), "== b")
print()
print("the inverse route gives the same answer here:")
print(" inv(A) @ b :", np.round(np.linalg.inv(A) @ b, 6))
print()
print("but solve is faster and numerically better - see the next step.")
Output
Why inverting is the wrong habit
It costs more and loses precision, and on a near-singular matrix it loses a lot.
example_04.pyNumPy
import numpy as np
# A Hilbert matrix: famously ill-conditioned, and we know the true answer.
n = 12
H = 1.0 / (np.arange(n)[:, None] + np.arange(n)[None, :] + 1)
x_true = np.ones(n)
b = H @ x_true
x1 = np.linalg.solve(H, b)
x2 = np.linalg.inv(H) @ b
e1 = float(np.abs(x1 - x_true).max())
e2 = float(np.abs(x2 - x_true).max())
print("condition number : %.2e <- badly conditioned" % np.linalg.cond(H))
print()
print("error with solve : %.3e" % e1)
print("error with inv : %.3e" % e2)
print("inv is %.0fx worse" % (e2 / e1))
print()
print("Both answers should be all ones. Neither is, because the problem")
print("itself is unstable - but inv throws away noticeably more.")
print("Rule: if you are about to write inv(A) @ b, write solve(A, b).")
Output
Determinant, rank and singular matrices
How to tell whether a system has a unique solution before you trust one.
example_05.pyNumPy
import numpy as np
good = np.array([[2.0, 1.0], [1.0, 3.0]])
bad = np.array([[1.0, 2.0], [2.0, 4.0]]) # row 2 = 2 x row 1
for name, M in [("good", good), ("singular", bad)]:
print("%-9s det=%8.4f rank=%d" % (name, np.linalg.det(M),
np.linalg.matrix_rank(M)))
print()
try:
np.linalg.solve(bad, np.array([1.0, 2.0]))
except np.linalg.LinAlgError as e:
print("solve(singular) ->", type(e).__name__, ":", e)
print()
print("But a determinant test is unreliable. This matrix is perfectly")
print("well behaved, and its determinant is nearly zero purely because")
print("the entries are small:")
tiny = 0.001 * np.eye(4)
print(" det = %.1e <- looks singular" % np.linalg.det(tiny))
print(" rank = %d of 4 <- it is not" % np.linalg.matrix_rank(tiny))
print(" cond = %.2f <- and it is beautifully conditioned"
% np.linalg.cond(tiny))
print()
print("Use matrix_rank to ask 'solvable?' and cond to ask 'trustworthy?'")
Output
Norms, eigenvalues and least squares
The three other things people actually reach for.
example_06.pyNumPy
import numpy as np
v = np.array([3.0, 4.0])
print("norm (length) :", float(np.linalg.norm(v)))
print("distance a to b :", float(np.linalg.norm(np.array([1.0, 1.0]) - v)))
S = np.array([[2.0, 1.0],
[1.0, 2.0]])
print()
print("symmetric matrix -> eigvalsh, ascending and real:")
print(" eigenvalues:", np.round(np.linalg.eigvalsh(S), 6))
print()
print("least squares fit of y = m x + c:")
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.1, 4.9, 7.0])
M = np.column_stack([x, np.ones_like(x)])
sol = np.linalg.lstsq(M, y, rcond=None)[0]
print(" slope %.4f intercept %.4f" % (sol[0], sol[1]))
Output
Shapes
(m, k) @ (k, n) gives (m, n). The inner dims must match and they cancel.
A 1-D array is special-cased: it's promoted to whichever orientation makes the product valid. Then the added dim is removed from the result. (2,3) @ (3,) gives (2,), treating the vector as a column. (3,) @ (3,4) gives (4,), treating it as a row.
That convenience means you rarely need to reshape vectors. It also means a shape error can come from a vector being the wrong length rather than the wrong orientation. The error message names both shapes. This Is often enough.
Solving systems
To solve Ax = b, use np.linalg.solve(A, b).
The obvious option, np.linalg.inv(A) @ b, gives the same answer on well-behaved problems and is wrong as a habit for two reasons.
Cost. finding an inverse is roughly the work of solving the system n times. solve factors the mat once and does one back-substitution.
score. The inverse adds rounding error. Then timesing by it adds more. solve avoids forming the inverse at all. On a well-conditioned mat the gap is small. On an ill-conditioned one it's the gap between a usable answer and noise.
The rule is simple enough to apply without thinking: if you're about to write inv(A) @ b, write solve(A, b).
Genuine uses for inv exist. You need the inverse itself, for a co-var mat or an analytical derivation — but they're rarer than the code you'll see suggests.
Is it solvable?
A singular mat has no unique answer, and solve raises LinAlgError.
The instinct is to check det(A) == 0 first. Don't. Determinants of floating-point mats are almost never just zero. The determinant scales badly. A well-conditioned mat can have a tiny determinant simply because it's large.
np.linalg.matrix_rank(A) answers the question right: full rank means solvable.
np.linalg.cond(A) measures how much the answer amplifies input error. A rule number near 1 is great. Anything above roughly 1e10 means you shouldn't trust the low-order digits of the result, whichever method made it.
Norms
np.linalg.norm(v) gives the Euclidean length of a vector. norm(a - b) the distance between two points.
It takes an axis, so norm(X, axis=1) gives the length of each row at once. The usual way to find distances for a whole dataset without looping.
Other norms are on hand via ord: ord=1 for the sum of abs values, ord=np.inf for the max.
Eigenvalues
np.linalg.eig works on any square mat and returns complex results in unspecified order.
For a even mat — which covers most cases that arise since co-var and Gram mats are even — use eigh or eigvalsh. They exploit the symmetry, are faster and more accurate. Return real eigen-vals in ascending order.
That determinism matters: code that indexes into the output of eig and assumes an order is leaning on a thing not sure.
Least squares
np.linalg.lstsq(M, y, rcond=None) fits an overdetermined system. More equations than unknowns. By cutting squared error.
Building the design mat with column_stack is the whole setup. For a straight line, the columns are x and a column of ones. The answer is slope and intercept.
Pass rcond=None by name to get the current action and silence the future-warning about the old default.
Batched matrix multiplication
matmul treats the last two axes as the mat and broadcasts all before them.
That makes (batch, n, k) @ (batch, k, m) a batch of mat products, found in one call with no loop. (batch, n, k) @ (k, m) also works, applying the same mat to each item in the batch.
This is why the batch axis goes first by rule. It puts the mat axes where matmul expects them. The whole thing composes without transposes.
np.dot doesn't do this. On arrays with more than two dims it follows an older rule that is rarely what anyone wants. @ and matmul are the ones to use.
Most of np.linalg also broadcasts over leading axes: solve, inv, det and the breakdowns all accept a stack of mats and return a stack of results.
solve versus lstsq
solve needs a square mat and gives the exact answer when one exists. It raises on a singular mat.
lstsq handles rectangular systems and finds the answer cutting squared error. It also handles rank-deficient systems by returning the min-norm answer rather than raising.
For an overdetermined system — more equations than unknowns. This Is what fitting a model to data looks like — lstsq is the right tool.
The tempting option is the normal equations: solve(A.T @ A, A.T @ b). It gives the same answer in exact math and is by number worse. Forming A.T @ A squares the rule number. A problem that was marginally conditioned turns badly conditioned. lstsq uses an SVD-based approach that avoids this.
SVD, and what it is for
np.linalg.svd factors any mat into U @ diag(s) @ Vt.
The singular values s are the handy part for most hands-on purposes.
Rank. The number of singular values usefully above zero. This is what matrix_rank finds inside. It's why rank is a more solid singularity test than a determinant.
Conditioning. The ratio of largest to smallest singular value is the rule number.
Dimensionality reduction. Keeping the largest k singular values and their vectors gives the best rank-k guess of the mat. That is what PCA is. Truncated SVD is how it's found.
Pseudo-inverse.np.linalg.pinv uses SVD to invert what can be inverted and skip what can't. This Is how least squares handles rank deficiency.
svd(a, full_matrices=False) returns the economy version. This Is smaller and often what you want for a tall mat.
Other decompositions
np.linalg.qr factors into an perp and an upper-tri mat. It's the by number sound way to solve least squares by hand. It's what many methods use inside.
np.linalg.cholesky factors a even plus-definite mat into a lower tri factor times its transpose. It's roughly twice as fast as a general split and is the standard tool for co-var mats. Sampling from a multivariate normal, or checking a Gaussian likelihood.
Cholesky raises if the mat isn't positive definite. This Is a handy diagnostic in itself. A co-var mat that fails Cholesky is telling you a thing is wrong with how it was guessed.
Conditioning in practice
The rule number is the amplification factor from input error to output error.
A rule number of 10 means roughly one decimal digit lost. 1e8 means about eight. This On float64's sixteen leaves eight. 1e16 means all of them. The answer is noise even so of the method.
Two things follow.
**Check cond before trusting a answer** from a mat built out of measured data. Nothing warns you.
Ill-conditioning is a property of the problem, not the code. No solver rescues a rule number of 1e16. The fix is upstream. Rescale vars so their magnitudes are close, drop collinear columns, or add fine. Which is just what ridge regression does, and why it works.
Scaling is the most often overlooked of those. A design mat with one column in metres and another in nanometres is badly conditioned for no reason other than units.
What to reach for
**@** for products, including batched.
**solve** for square systems, never inv.
**lstsq** for fitting and overdetermined systems, never the normal equations.
**eigvalsh / eigh** for even mats. eig returns complex values in unspecified order.
**cholesky** for co-var mats.
**svd** when you need rank, conditioning, or a low-rank guess.
**cond and matrix_rank** before believing any of it.
And for anything beyond this. Sparse mats, iterative solvers, specialised breakdowns — SciPy's scipy.linalg and scipy.sparse extend the same interface and are where the rest of the toolkit lives.
Common mistakes with Linear Algebra
**Using * where @ was meant.** Both are legal on square mats and give other answers of the same shape. Nothing warns. This is the most frequent linear algebra bug in NumPy code.
**Writing inv(A) @ b.** Slower and less accurate than solve(A, b), for no gain.
Using the normal equations for least squares.solve(A.T @ A, A.T @ b) squares the rule number. lstsq doesn't.
**Testing det(A) == 0.** Floating-point determinants are almost never just zero. The determinant scales with the mat. matrix_rank answers the question right.
**Indexing into eig output given an order.** There's none. eigh and eigvalsh return ascending real values for even mats. This Is what most code really needs.
**Expecting .T to make a column vector.** On a 1-D array it does nothing. v[:, None].
**Using np.dot on more than two dims.** It follows an older rule that is rarely what anyone wants. @ or matmul broadcast the leading axes properly.
Reading the shapes in an error
matmul errors name the "core dims". This Are the last two axes.
"Input operand 1 has a mismatch in its core dim 0" means the first of the last two axes of the second operand doesn't match what the first operand's last axis needs. The k in (m,k) @ (k,n).
The diagnosis is to print both shapes and check which axis is supposed to be shared. Nearly always one of the two arrays is transposed scaled to the intent. The fix is a .T on the correct side rather than a reshape.
For a 1-D operand, note it's promoted to whichever orientation makes the product valid. A shape error involving one often means the wrong length rather than the wrong orientation.
Questions people ask about Linear Algebra
What is the difference between * and @? The star multiplies element by element. The @ operator is matrix multiplication, which is the one you want for linear algebra.
Should I use np.linalg.inv to solve Ax = b? No. Use np.linalg.solve(A, b), which is faster and more accurate than forming the inverse and multiplying.
How do I check whether a matrix is singular? Look at np.linalg.matrix_rank or the condition number from np.linalg.cond. A determinant near zero is a weaker test because it depends on scale.
What do eig and eigh differ in? eigh is for symmetric or Hermitian matrices. It is faster, returns real eigenvalues in ascending order, and is the right choice for covariance matrices.
How do I fit a least-squares line? Use np.linalg.lstsq(A, y, rcond=None). It handles non-square and rank-deficient systems that solve cannot.
What does np.dot do compared with @? For 1D and 2D arrays they agree. For higher dimensions they differ, and matmul (@) broadcasts over leading axes, which is usually what you want.
Where SciPy takes over
np.linalg covers the dense, general cases well. a few things sit just outside it.
Sparse mats.scipy.sparse stores only nonzeros and has its own solvers. For a mat that is mostly zeros. Graphs, finite-element meshes, term-doc mats. The gap isn't a constant factor but the gap between fitting in memory and not.
Iterative solvers. For very large systems, scipy.sparse.linalg offers methods that rough a answer without factoring the mat.
More breakdowns.scipy.linalg has LU, Schur, banded and tri solvers. Often more options on the ones NumPy also gives.
mat functions.expm for the mat steep and friends live in SciPy.
The interfaces are on purpose similar. Moving across is often a change of import rather than a change of approach.
Check yourself
0 of 4
Answer without scrolling back up.
What is the difference between `a * b` and `a @ b` for two 2x2 arrays?
Both are legal and give different answers of the same shape, with no warning. It is the most common linear algebra bug in NumPy code.
Why prefer `np.linalg.solve(A, b)` over `np.linalg.inv(A) @ b`?
solve factors once and back-substitutes. On an ill-conditioned matrix the accuracy gap is the difference between an answer and noise.
How should you check whether a matrix is singular?
Floating-point determinants are almost never exactly zero, and a determinant can be tiny purely because the matrix is large.
Why use `eigvalsh` rather than `eig` for a symmetric matrix?
`eig` returns complex values in unspecified order, so any code indexing its output is relying on something not guaranteed.
Cheat sheet
Linear Algebra
Multiply matrices, solve linear systems and compute eigenvalues with NumPy. See why you should call solve instead of computing the inverse.
Ashish Jangra builds and maintains VizLearn. Every module here is written and the visualisation behind it hand-built, so the numbers in a readout come from the same code that draws the picture. Corrections are genuinely welcome and get priority over everything else — if a page states something wrong, or an animation misrepresents what the algorithm does, get in touch.