scipy.linalg — more stable than numpy.linalg, lu, qr, schur: Visual walkthrough

Step 1 — What "least squares" even asks

WHAT. We have more equations than unknowns. Picture a machine (a matrix $A$ ) that eats a small arrow $x$ and spits out an arrow in a big space. We have a target arrow $b$ out there. We ask: which input $x$ makes the output $A x$ land as close to $b$ as possible?

WHY this question. In real data (fitting a line to 100 noisy points) there is no exact solution — the points don't lie on one line. So "solve $A x = b$ exactly" is impossible; the honest question becomes "get as close as you can."

PICTURE. Below, $b$ is the orange target. As $x$ varies, $A x$ can only reach points on the blue column space — the flat plane of everything the machine can produce. The closest reachable point is $A x^{⋆}$ (green). The leftover gap is the red residual.

Figure — scipy.linalg — more stable than numpy.linalg, lu, qr, schur

So the whole problem is: $x min ∥ b - A x ∥$ — make the red arrow as short as possible.

Step 2 — The geometric answer: drop a perpendicular

WHAT. Before any algebra, here is the answer in one idea: the shortest red arrow from the target $b$ down to the blue plane is the perpendicular one.

WHY. Think of standing at point $b$ above a flat floor (the column space). The closest floor point to you is straight down — any slanted path is longer. That "straight down" is the perpendicular drop. So at the best $x^{⋆}$ , the residual $r$ must be perpendicular to the whole plane.

PICTURE. The green foot of the perpendicular is $A x^{⋆}$ ; the red $r$ meets the plane at a right angle (see the little square).

Everything after this is just: how do we compute that perpendicular drop cleanly?

Step 3 — Length is what we must protect (enter orthogonal matrices)

WHAT. Our whole problem is about length ( $∥ b - A x ∥$ ). So the perfect tool is a transformation that does not change lengths. That tool is an orthogonal matrix $Q$ .

WHY this tool and not another? We are free to look at the problem "through a rotated pair of glasses." But we must not distort the picture — if the glasses stretched things, the shortest arrow through them would be a different arrow, and we'd solve the wrong problem. Orthogonal matrices are exactly the rotations/reflections: they spin and flip space but never stretch it.

PICTURE. Left: an arrow $v$ and its length. Right: after applying $Q$ (a rotation), the arrow points elsewhere but the ruler reads the same length.

Step 4 — Build $Q$ and $R$ : straightening the columns of $A$

WHAT. We factor the machine as $A = QR$ :

$Q$ = an orthonormal set of directions (clean perpendicular axes),
$R$ = an upper-triangular recipe saying how each original column of $A$ is mixed from those axes.

WHY. $A$ 's columns are usually slanted and unequal in length — awkward. Gram–Schmidt straightens them into perpendicular unit arrows one at a time, and bookkeeps the straightening in $R$ . We separate the machine into "a clean rotation" ( $Q$ ) times "a simple triangular stretch" ( $R$ ).

PICTURE. Two slanted columns $a_{1}, a_{2}$ (blue) become perpendicular unit vectors $q_{1}, q_{2}$ (green). $R$ (orange) is the small triangular table of coefficients: $a_{1} = r_{11} q_{1}$ , and $a_{2} = r_{12} q_{1} + r_{22} q_{2}$ .

Step 5 — Look at the problem through $Q$ 's glasses

WHAT. Substitute $A = QR$ into the residual, then rotate the whole picture by $Q^{⊤}$ (the reverse rotation).

WHY. From Step 3, $Q^{⊤}$ is orthogonal too, so it doesn't change any lengths — the shortest arrow stays the shortest arrow. But it aligns everything to $Q$ 's clean axes, where the algebra falls apart into something trivial.

PICTURE. Same target/plane as before, but now viewed after applying $Q^{⊤}$ : the slanted column space becomes axis-aligned, and the residual measurement is unchanged (same red length).

Step 6 — Read off the answer: $R x = Q^{⊤} b$

WHAT. In the rotated view, part of $Q^{⊤} b$ lands inside the reachable space (the top $n$ coordinates, matched by $R x$ ) and part lands outside it (the leftover bottom coordinates, which $R x$ can never touch).

WHY. We can drive the reachable part exactly to zero by choosing $x$ so that $R x$ equals the top part of $Q^{⊤} b$ . The unreachable bottom part is the irreducible error — the length of the perpendicular from Step 2. Nothing we do to $x$ shrinks it.

PICTURE. $Q^{⊤} b$ split into a "reachable" top block (green, cancellable) and an "unreachable" bottom block (red, the minimum residual). The best $x$ zeroes the green.

Step 7 — Edge and degenerate cases (never leave a scenario unshown)

WHAT & WHY. A derivation you can trust must survive the weird inputs.

The figure shows all three: (A) perpendicular meets a point → exact; (B) a collapsed axis → non-unique foot; (C) target on the plane → zero residual.

The one-picture summary

The full journey on one canvas: target $b$ (orange) → drop the perpendicular to the column space (green foot, red residual) → rotate by $Q^{⊤}$ so the space becomes axis-aligned → the reachable part is cancelled by $R x$ , leaving the unreachable block as the minimum miss → solve the triangular $R x = Q^{⊤} b$ by back-substitution.

Recall Feynman retelling — say it to a friend

You've got a target dot floating above a flat table. You can only place your pencil on the table, and you want to be as close to the target as possible. Obviously: stand directly beneath it — straight down. That "straight down" is the least-squares answer, and the leftover height is the error you simply can't remove.

Doing "straight down" with messy slanted table-legs is annoying. So we first straighten the table's edges into clean perpendicular rulers (that's $Q$ ), writing down how the old edges were made from the new ones (that's $R$ ). Because rotating the whole scene with clean rulers never changes any distances, the closest-point problem stays identical but now reads off instantly: match the reachable part exactly, accept the unreachable leftover. What's left is a tiny triangular puzzle $R x = Q^{⊤} b$ you solve by walking one row at a time.

The tempting shortcut — squaring $A$ into $A^{⊤} A$ — secretly doubles your error, because it squares the condition number. QR refuses to do that. Same answer on paper, far more trustworthy on a computer. In practice you never hand-crank it: x, res, rank, sv = scipy.linalg.lstsq(A, b).

Prereqs revisited: Least Squares, Gram-Schmidt Orthogonalization, Condition Number, LAPACK and BLAS, numpy.linalg.