Numpy linalg.lstsq – Return the least-squares solution to a linear matrix equation

I had three points and a line to fit, and I kept writing normal equations by hand. I expected residuals to always be a number, so the empty array on my first exact system felt like a bug until I read the four returns.

When a linear system has no exact solution

A system Ax = b has an exact solution only when b lies in the column space of A. When you have more equations than unknowns, that rarely happens, so the best you can do is minimize the squared error, which is the least-squares problem. I ran a tiny 3 by 2 system and the solver gave me the line that missed each point by the smallest total squared distance.

Least squares chooses x that minimizes sum((Ax – b)**2). That choice is also what numpy.linalg.lstsq computes via SVD, so you get not just x but the diagnostics that tell you whether to trust it.

import numpy as np

A = np.array([[1,1],[1,2],[1,3]], dtype=float)
b = np.array([1,2,2], dtype=float)
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print(x)          # [0.66666667 0.5]
print(residuals)  # [0.166...]
print(rank, s)    # 2 [4.07 0.6]

What lstsq needs before you call it

You need NumPy and the right shapes, plus a conscious choice for rcond. I tested on NumPy 2.4.6 with Python 3.11.16, so rcond=None uses the current default that replaces the old warn behavior. A is m by n, b is m or m by k.

import numpy as np
print(np.__version__)  # 2.4.6 on my run

# rcond=None is the current safe default
A = np.array([[1,1],[1,1.0001]])
b = np.array([2,2.0001])
x, *_ = np.linalg.lstsq(A, b, rcond=None)
print(x)

How to solve with lstsq and read the four answers

lstsq returns four things, and each answers a distinct question, so unpack them by name.

Solve the basic overdetermined case

The regression-style case has more rows than columns. The solution x is the best fit, and residuals tells you how much error remains.

import numpy as np

A = np.array([[1,1],[1,2],[1,3]], dtype=float)
b = np.array([1,2,2], dtype=float)
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print(f"x={x}")
print(f"Ax={A @ x}")
print(f"residuals={residuals}")  # sum of squared errors
print(f"rank={rank}, s={s}")

I assumed residuals would always contain a number. Instead the overdetermined case gave 0.166…, which is the actual fit quality I now check before trusting x.

Handle exact and underdetermined systems

When the system is exactly determined or underdetermined, residuals comes back empty. That is not an error, it is the API telling you there was no residual to compute.

import numpy as np

# Exact 2x2
A2 = np.array([[2,1],[1,3]], dtype=float)
b2 = np.array([8,13], dtype=float)
x2, res2, rank2, s2 = np.linalg.lstsq(A2, b2, rcond=None)
print(x2, res2)  # [2.2 3.6] [] empty because exact

# Underdetermined 1x3
A3 = np.array([[1,1,1]], dtype=float)
b3 = np.array([3], dtype=float)
x3, res3, rank3, s3 = np.linalg.lstsq(A3, b3, rcond=None)
print(x3, res3, rank3)  # [1. 1. 1.] [] 1

The empty residuals caught me once. I now branch on len(residuals) rather than treating it as zero, because zero and empty mean different things here.

Run a real regression with lstsq

For y = slope*x + intercept, build A as [x, ones] and solve. This is where lstsq replaces hand-rolled formulas.

import numpy as np

rng = np.random.default_rng(0)
X = np.linspace(0,10,20)
y = 2*X + 1 + rng.normal(0,0.5,20)
A = np.column_stack([X, np.ones_like(X)])
coeff, residuals, rank, s = np.linalg.lstsq(A, y, rcond=None)
print(coeff)  # [1.97 1.04] close to true [2 1]

# compare with pinv, same result via different path
print(np.allclose(coeff, np.linalg.pinv(A) @ y))  # True

The fitted slope came back 1.97 against a true 2, so with noise the error was visible but bounded, which is exactly what rank 2 and the singular values promised.

Where lstsq answers need checking

Three checks keep you from misreading the output.

import numpy as np

# Rank deficient -> small s signals trouble
A = np.array([[1,1],[1,1.0001]])
print(np.linalg.lstsq(A, [2,2.0001], rcond=None)[3])  # s around [2.0, 5e-05]

# rcond filters tiny singular values
x_default, *_ = np.linalg.lstsq(A, [2,2.0001], rcond=None)
x_strict, *_ = np.linalg.lstsq(A, [2,2.0001], rcond=1e-3)
print(x_default, x_strict)

# Always check residuals emptiness
A = np.array([[1,1],[1,2],[1,3]], float)
b = np.array([1,2,2], float)
_, res, _, _ = np.linalg.lstsq(A,b,rcond=None)
print("has residuals:", len(res)>0, res)

• Empty residuals means exact or underdetermined, not zero error
• Tiny s values mean A is nearly rank deficient, so x is sensitive
• Set rcond only when you know the noise floor, otherwise keep None

What you now have

You have a runnable pattern for overdetermined, exact, and underdetermined cases plus the rule for the four returns. I keep the unpack x, residuals, rank, s and check len(residuals) and s before I trust x, because those two diagnostics are where the silent mistakes live.

# minimal template I now copy
import numpy as np
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
ok = len(residuals)>0 or rank==A.shape[1]
print(rank, s[:2])

FAQ