Don't Invert That Matrix

In data analysis you are lucky if your problem can be boiled down to solving a linear system of equations $Ax = b$. This comes up routinely in spectroscopy, particularly for nuisances or systematics: fitting a continuum to a spectrum, fitting a low-order polynomial to some calibration data, solving for the coefficients of a linear model.

In these situations, your first instinct might be to reach for the inverse:

x = np.linalg.inv(A) @ b

It reads exactly like the maths, $x = A^{-1}b$. It gives you an answer. For the simplest and well-posed problems, it works. So you keep doing that while you expand to work on more complex problems. Then you find it still gives you an answer, but that answer is wrong. Only then do you learn that explicitly computing an inverse like this is almost never the right way to solve for $x$.

The naive move: normal equations

Let's say you have some data and model parameters. We'll assume you have more data than parameters. $A$ is a matrix with shape $m \times n$ with $m > n$, so there are a couple of problems. The matrix is not square. And there's no exact solution to $Ax = b$, only a least-squares one. The textbook derivation says the minimizer of $\|Ax - b\|^2$ satisfies the so-called normal equations

and $A^T A$ is square, so now you can invert it:

x = np.linalg.inv(A.T @ A) @ (A.T @ b)

This is still the thing not to do. Forming $A^T A$ squares the condition number of the problem,

so if $A$ is even mildly ill-conditioned, $A^T A$ is much worse, and you're inverting a matrix that is numerically far closer to singular than the one you started with.

To remind you about condition numbers: a matrix doesn't just sit there holding numbers, it does something -- you feed it a vector and it stretches it. It stretches some directions a lot and others barely at all, and the singular values are exactly those stretch factors. The condition number is the ratio of the biggest stretch to the smallest one. (It is tempting to eyeball the biggest and smallest entries of the matrix instead, but that will mislead you: a matrix whose entries are all close to 1 can still have a condition number of $10^{15}$, if two of its rows are nearly the same as each other.)

Why care? Because solving $Ax = b$ means running that stretching backwards. If the matrix squashed some direction down to almost nothing on the way in, then undoing it on the way out means blowing that direction back up by an enormous factor -- and whatever rounding error happened to be sitting in it gets blown up right along with it. That is the "wiggle room" you are spending. A double gives you about 16 digits; a condition number of $10^{12}$ has spent twelve of them before you have done anything at all, and you are left doing your arithmetic with the four that survive. When your condition number is large, you are more susceptible to having errors in your results, or hitting numerical overflow.

Here's a concrete case: fitting a degree-10 polynomial (an 11-parameter model) to 60 points with numpy.vander, which is exactly the kind of design matrix you get from polynomial continuum fitting.

import numpy as np
from scipy.linalg import lstsq

rng = np.random.default_rng(0)
x = np.linspace(0, 1, 60)
A = np.vander(x, 11, increasing=True)

true_coeffs = rng.normal(size=11)
b = A @ true_coeffs  # noiseless, so any error is pure rounding error

x_normal = np.linalg.inv(A.T @ A) @ (A.T @ b)
x_lstsq, *_ = lstsq(A, b)

The condition numbers already tell the story:

cond(A)       = 2.05e+07
cond(A.T @ A) = 4.21e+14

Big condition number bad, small condition number good. Someone told me it's like something called 'golf', if that helps you. Since we know the true values true_coeffs, we can compare the coefficients we recovered against the true values:

Remember: these data are noiseless. There is an exact answer, and lstsq finds it to near machine precision. The normal-equations approach, on the same data, is already wrong in the second decimal place, purely from rounding error introduced by squaring the condition number. lstsq avoids this because it works directly on $A$ — using methods like singular value decomposition (SVD) or QR factorization — which means it never needs to form the product $A^T A$.

Even for square systems, skip the inverse

Even when $A$ is square and least squares doesn't enter into it, inv(A) @ b is still not the move you want.¹

numpy.linalg.solve factorises $A$ (through LU decomposition, then forward/back substitution) and never computes the other $n-1$ columns of $A^{-1}$ that you were going to throw away. That means it is also faster:

You can see np.linalg.solve is around four times faster for the same answer. With a naive inv(M) you factorise the whole matrix to get $n^2$ numbers back when you only needed $n$. Speed is the less important half, though. Explicit inversion also amplifies error. Let's build a deliberately ill-conditioned $200 \times 200$ system:

rng = np.random.default_rng(2)

n = 200
U, _, Vt = np.linalg.svd(rng.standard_normal((n, n)))
s = np.logspace(0, -12, n)      # kappa = s[0] / s[-1] = 1e12
A = (U * s) @ Vt

x_true = rng.standard_normal(n)
b = A @ x_true

x_inv = np.linalg.inv(A) @ b
x_solve = np.linalg.solve(A, b)

The condition number will be about $\kappa(A) \approx 10^{12}$ — not exotic, this happens routinely with near-degenerate design matrices, badly scaled features, or nearly-collinear basis functions — so now let's check how it compares against a known true solution:

cond(A) = 1.0e+12
inv  : relative error 1.8e-04   residual 2.3e-05
solve: relative error 9.3e-06   residual 1.2e-15

Note the two columns disagree about how big the win is. Note the two columns disagree about how big the win is! solve is giving us a model that is about ten orders of magnitude better in terms of residuals (e.g., $\chi^2$) than what we recover with inv. When we are doing a direct inv we are putting the computational emphasis on "knowing $A^{-1}$ well" instead of putting the emphasis where we should be: "knowing $x$ well".

That gap is the condition number doing its work: at $\kappa \approx 10^{12}$ a tiny residual no longer guarantees an accurate $x$. solve isn't immune to ill-conditioning — nothing is — it just doesn't add insult to injury by materialising the inverse first.

When you have constraints: lsq_linear

lstsq and solve will happily hand you negative coefficients, or coefficients outside whatever physical range makes sense: a negative flux, a negative abundance, an unphysical column density. If your parameters are genuinely constrained, scipy.optimize.lsq_linear solves the same linear least-squares problem subject to bounds, using a trust-region-reflective algorithm, and still never forms an inverse:

from scipy.optimize import lsq_linear

res = lsq_linear(A, b, bounds=(0, np.inf))
res.x        # all >= 0
res.success  # True

On a random $2000 \times 50$ problem, the unconstrained fit puts 28 of the 50 coefficients below zero; bounding them at zero leaves none negative and moves the residual from 44.24 to 44.54. That is the trade you are making, and it's worth saying out loud: imposing bounds always makes the residual worse, because you're deliberately ruling out the unconstrained optimum. lsq_linear is for when the constraint is real physics, not for tidying up an ugly-looking fit.

When the model itself is non-linear: least_squares

lsq_linear is still solving a linear problem — the model is linear in its parameters, just bounded. If your model isn't linear in its parameters (a Voigt profile, an exponential decay, anything you'd otherwise point curve_fit at) you need scipy.optimize.least_squares. You supply a residual function, and ideally its Jacobian. It supports bounds too, and unlike the others it can do robust loss functions (loss="soft_l1", "huber") to downweight outliers like cosmic-ray hits.

But don't use it by default: on that same $2000 \times 50$ linear problem, all three agree on the answer to $8\times10^{-17}$, but not on the cost of getting there:

Two things to take from that. First, least_squares is doing iterative nonlinear optimisation on a problem with a closed-form solution — it has no way to know your residual is linear in $x$, so it pays for generality you don't need. Second, look at the last two rows: most of the penalty isn't the iteration, it's the finite-difference Jacobian. If you are solving something genuinely nonlinear, writing the Jacobian by hand is the single highest-leverage thing you can do.

There's also a correctness argument, not just a speed one. A linear least-squares problem has exactly one global optimum. A nonlinear solver started somewhere unhelpful does not come with that guarantee.

When you actually do want the inverse

There are legitimate cases, and the honest one is uncertainties: $(J^T J)^{-1}$ is a covariance matrix, and you want its entries, not a solution vector. Even then, don't call inv on the normal equations. Since $J^T J$ is a covariance matrix, then you know that it is symmetric and positive-definite (or non-negative-definite), so you should use a Cholesky factorisation (scipy.linalg.cho_factor and cho_solve) and solve for the columns you need — or just the diagonal if you only need the per-parameter error bars. It's the same principle as everywhere else in this post: use the factorisation that matches the structure of your problem instead of the most general and most expensive tool available.

The decision

If you find yourself writing inv(...) @ ..., you are almost always one step away from a more direct — and more numerically honest — way to get the same $x$.

inv(A) @ b, and reaching for least_squares on a problem that's actually linear, are the same mistake wearing different clothes: using the most general and most expensive tool for a problem that has a cheaper, more specific solution.

Reproduce this

Every number above comes from the script below. Each section seeds its own generator, so you can run or edit them in any order and get reproducible results:

import time

import numpy as np
from scipy.linalg import lstsq
from scipy.optimize import least_squares, lsq_linear


def best_of(f, n=7):
    """Best-of-n wall clock, in milliseconds."""
    times = []
    for _ in range(n):
        t0 = time.perf_counter()
        f()
        times.append(time.perf_counter() - t0)
    return min(times) * 1e3


def section(title):
    print(f"\n{title}\n{'-' * len(title)}")


# ---------------------------------------------------------------------------
section("1. Normal equations vs lstsq (degree-10 polynomial, 60 points)")

rng = np.random.default_rng(0)

x = np.linspace(0, 1, 60)
A = np.vander(x, 11, increasing=True)

true_coeffs = rng.normal(size=11)
b = A @ true_coeffs  # noiseless, so any error is pure rounding error

x_normal = np.linalg.inv(A.T @ A) @ (A.T @ b)
x_lstsq, *_ = lstsq(A, b)

print(f"cond(A)       = {np.linalg.cond(A):.2e}")
print(f"cond(A.T @ A) = {np.linalg.cond(A.T @ A):.2e}")
print(f"normal eqns : coefficient error {np.linalg.norm(x_normal - true_coeffs):.1e}"
      f"   residual {np.linalg.norm(A @ x_normal - b):.1e}")
print(f"lstsq       : coefficient error {np.linalg.norm(x_lstsq - true_coeffs):.1e}"
      f"   residual {np.linalg.norm(A @ x_lstsq - b):.1e}")

# ---------------------------------------------------------------------------
section("2. inv vs solve, speed (square, well-conditioned)")

rng = np.random.default_rng(1)

for n in (500, 1500):
    M = rng.normal(size=(n, n))
    v = rng.normal(size=n)
    t_inv = best_of(lambda: np.linalg.inv(M) @ v)
    t_solve = best_of(lambda: np.linalg.solve(M, v))
    print(f"n = {n:5d}:  inv(M) @ v {t_inv:7.1f} ms"
          f"   solve(M, v) {t_solve:6.1f} ms   ({t_inv / t_solve:.1f}x)")

# ---------------------------------------------------------------------------
section("3. inv vs solve, accuracy (ill-conditioned, n = 200)")

# The matrix is *constructed* to have cond(A) = 1e12, not found in the wild:
# take a random matrix, then replace its singular values with a geometric
# sweep spanning twelve decades. This is the honest way to hit a target
# condition number; s[0] * logspace(0, -12, n) gives kappa = s[0]/s[-1] = 1e12.
rng = np.random.default_rng(2)

n = 200
U, _, Vt = np.linalg.svd(rng.standard_normal((n, n)))
s = np.logspace(0, -12, n)
A = (U * s) @ Vt

x_true = rng.standard_normal(n)
b = A @ x_true

x_inv = np.linalg.inv(A) @ b
x_solve = np.linalg.solve(A, b)

rel = lambda xh: np.linalg.norm(xh - x_true) / np.linalg.norm(x_true)
print(f"cond(A) = {np.linalg.cond(A):.1e}")
print(f"inv   : relative error {rel(x_inv):.1e}"
      f"   residual {np.linalg.norm(A @ x_inv - b):.1e}")
print(f"solve : relative error {rel(x_solve):.1e}"
      f"   residual {np.linalg.norm(A @ x_solve - b):.1e}")

# ---------------------------------------------------------------------------
section("4. Bounds bind, and cost you residual (2000 x 50)")

rng = np.random.default_rng(3)

A = rng.standard_normal((2000, 50))
b = rng.standard_normal(2000)

x_free, *_ = lstsq(A, b)
res_bounded = lsq_linear(A, b, bounds=(0, np.inf))

print(f"unconstrained : {(x_free < 0).sum()} of 50 coefficients negative,"
      f" residual {np.linalg.norm(A @ x_free - b):.2f}")
print(f"bounded >= 0  : {(res_bounded.x < -1e-12).sum()} negative,"
      f" residual {np.linalg.norm(A @ res_bounded.x - b):.2f}")

# ---------------------------------------------------------------------------
section("5. lstsq vs lsq_linear vs least_squares on the same linear problem")

# Same A, b as section 4 — deliberately, so the timings and the bounds result
# describe one problem rather than two.
x0 = np.zeros(50)
residual_fn = lambda p: A @ p - b
jacobian_fn = lambda p: A

t_lstsq = best_of(lambda: lstsq(A, b))
t_lsq_linear = best_of(lambda: lsq_linear(A, b))
t_ls_jac = best_of(lambda: least_squares(residual_fn, x0, jac=jacobian_fn))
t_ls_fd = best_of(lambda: least_squares(residual_fn, x0), n=3)

print(f"lstsq                                {t_lstsq:6.1f} ms")
print(f"lsq_linear                           {t_lsq_linear:6.1f} ms"
      f"   ({t_lsq_linear / t_lstsq:.1f}x)")
print(f"least_squares, analytic Jacobian     {t_ls_jac:6.1f} ms"
      f"   ({t_ls_jac / t_lstsq:.1f}x)")
print(f"least_squares, finite differences    {t_ls_fd:6.1f} ms"
      f"   ({t_ls_fd / t_lstsq:.1f}x)")

x_ls = least_squares(residual_fn, x0, jac=jacobian_fn).x
print(f"max |lstsq - least_squares| = {np.abs(x_free - x_ls).max():.1e}"
      "  (same answer)")