least squares via QR factorization
The QR route is the standard, dependable way to solve least squares — the one a good library uses by default. The idea is to rewrite the tall matrix A as a product Q R, where Q has orthonormal columns (perpendicular unit vectors, so Q^T Q = I) and R is upper triangular. Because Q's columns are orthonormal, multiplying by Q^T preserves lengths exactly — it just rotates and reflects space without stretching it — so it cannot amplify error. That is the whole trick: do the hard work with length-preserving operations.
Here is the method in plain steps. Factor A = Q R (thin QR: Q is m-by-n, R is n-by-n upper triangular). Minimizing ||A x - b||_2 = ||Q R x - b||_2; multiply inside the norm by Q^T (which does not change the length): the problem becomes ||R x - Q^T b||_2 on the relevant components. So you compute the vector c = Q^T b, then solve the small triangular system R x = c by back substitution — fast and stable. Crucially, you NEVER form A^T A, so you never square the condition number. The accuracy you get tracks kappa(A), not kappa(A)^2.
QR factorizations are built by Householder reflections (the most common, backward-stable for dense problems), Givens rotations (handy for sparse or streaming updates), or Gram-Schmidt (intuitive but needs the 'modified' variant to be reliable). QR costs about twice the flops of the normal equations (~2 m n^2 vs ~m n^2 + n^3/3), but the extra cost buys numerical safety. For full-rank problems QR is the recommended default; for rank-deficient or near-singular ones, step up to the SVD.
To fit our line, factor A (with rows (1,1),(1,2),(1,3)) into Q (3-by-2, orthonormal columns) times R (2-by-2 upper triangular). Then c = Q^T b is a 2-vector, and back-substituting R x = c yields x = (2/3, 1)^T — the same answer as the normal equations, but obtained without ever squaring kappa(A).
Orthonormal Q preserves length, so QR solves least squares without squaring the conditioning.
QR is the recommended default for full-rank least squares; it costs about twice the normal equations but is far more accurate on ill-conditioned data. It does not, however, handle rank deficiency gracefully — for that, use the SVD.