6 Linear Systems
Most of applied linear algebra reduces to one question: given \(\mathbf{A}\) and \(\mathbf{b}\), find \(\mathbf{x}\) with
\[ \mathbf{A}\mathbf{x} = \mathbf{b} \]
Fitting a regression, interpolating a curve, balancing a chemical equation, and computing the stationary distribution of a Markov chain are all this, wearing different clothes.
This chapter answers three things: when a solution exists, how to find it, and — the part most courses skip and every practitioner eventually needs — what to do when the answer is technically correct but useless.
6.1 Systems of linear equations
A linear equation in the unknowns \(x_1, \dots, x_n\) has the form
\[ a_1x_1 + a_2x_2 + \cdots + a_nx_n = b \]
Linear means each unknown appears alone, to the first power, multiplied by a constant. No \(x^2\), no \(x_1x_2\), no \(\sin x\). A system is several such equations that must hold at once:
\[ \begin{aligned} 2x_1 + \phantom{3}x_2 &= 5 \\ \phantom{2}x_1 + 3x_2 &= 10 \end{aligned} \]
Stack the coefficients into a matrix and the right-hand sides into a vector, and the whole system becomes one matrix equation:
\[ \underbrace{\begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{\mathbf{x}} = \underbrace{\begin{bmatrix} 5 \\ 10 \end{bmatrix}}_{\mathbf{b}} \]
That compression is the point. Everything from Chapter 5 now applies.
6.2 The geometric picture
There are two ways to see a system, and they answer different questions.
6.2.1 The row picture
Each equation is a line (in 2D) or a plane (in 3D). A solution is a point lying on all of them at once — an intersection.
Two lines in a plane usually cross exactly once, which is why a square system usually has exactly one solution. The exceptions are the interesting part (Section 6.5).
6.2.2 The column picture
Recall from Section 5.7 that \(\mathbf{A}\mathbf{x}\) is a linear combination of \(\mathbf{A}\)’s columns. So solving \(\mathbf{A}\mathbf{x} = \mathbf{b}\) asks:
How much of each column do I need to build \(\mathbf{b}\)?
The column picture is the more useful of the two, because it makes solvability obvious: a solution exists exactly when \(\mathbf{b}\) lies in the span of the columns. In \(\mathbb{R}^2\) with two independent columns, that span is the whole plane, so every \(\mathbf{b}\) is reachable.
6.3 Gaussian elimination
The algorithm that solves systems by hand, and — with care about which row goes where — the one your computer actually uses.
The idea: three operations change the equations without changing the solution set.
- Swap two rows.
- Multiply a row by a non-zero constant.
- Add a multiple of one row to another.
Use them to eliminate unknowns until the last equation has just one, then work backwards.
Worked example.
\[ \begin{aligned} x_1 + \phantom{2}x_2 + \phantom{2}x_3 &= 6 \\ 2x_1 + \phantom{2}x_2 - \phantom{2}x_3 &= 1 \\ \phantom{2}x_1 - \phantom{2}x_2 + 2x_3 &= 5 \end{aligned} \]
Write the coefficients and right-hand side together as an augmented matrix:
\[ \left[\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 2 & 1 & -1 & 1 \\ 1 & -1 & 2 & 5 \end{array}\right] \]
Eliminate \(x_1\) from rows 2 and 3, using row 1 as the pivot row. Subtract \(2 \times\) row 1 from row 2, and \(1 \times\) row 1 from row 3:
\[ \left[\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & -1 & -3 & -11 \\ 0 & -2 & 1 & -1 \end{array}\right] \]
Now eliminate \(x_2\) from row 3, using row 2 as the pivot. Subtract \(2 \times\) row 2:
\[ \left[\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & -1 & -3 & -11 \\ 0 & 0 & 7 & 21 \end{array}\right] \]
The last row now reads \(7x_3 = 21\), so \(x_3 = 3\). Back-substitute upward:
\[ \begin{aligned} -x_2 - 3(3) &= -11 &&\Rightarrow x_2 = 2 \\ x_1 + 2 + 3 &= 6 &&\Rightarrow x_1 = 1 \end{aligned} \]
The solution is \(\mathbf{x} = (1, 2, 3)\).
[1] 1 2 3The cost is \(O(n^3)\), the same order as matrix multiplication, and the same reason: about \(n\) elimination steps, each touching about \(n^2\) entries.
6.4 Row echelon form
The staircase shape elimination produces has a name. A matrix is in row echelon form when:
- every all-zero row sits at the bottom, and
- each row’s leading non-zero entry (its pivot) is strictly to the right of the pivot in the row above.
\[ \begin{bmatrix} \boxed{1} & 1 & 1 \\ 0 & \boxed{-1} & -3 \\ 0 & 0 & \boxed{7} \end{bmatrix} \]
Keep going — scale every pivot to 1 and clear the entries above each pivot too — and you reach reduced row echelon form, where the solution can be read straight off without back-substitution. It is tidier by hand and slower by machine, so software does not bother.
The number of pivots is the rank of the matrix, and it is the single most informative number about a linear system. We give it a proper treatment in Section 7.8.
6.5 Existence and uniqueness
A square system has exactly one of three fates, and the row picture shows all three.
6.5.1 Exactly one solution
The lines cross once. Algebraically: \(n\) pivots for \(n\) unknowns, \(\det(\mathbf{A}) \neq 0\), \(\mathbf{A}\) invertible. This is Figure 6.1.
6.5.2 No solution
The equations contradict each other. Two parallel lines never meet.
Elimination gives itself away here: a row becomes \(\begin{bmatrix}0 & 0 & \cdots\end{bmatrix}\) on the left with something non-zero on the right — the equation \(0 = 1\). Such a system is inconsistent.
In the column picture, \(\mathbf{b}\) simply does not lie in the span of the columns.
6.5.3 Infinitely many solutions
The equations are redundant — \(x_1 + x_2 = 1\) and \(2x_1 + 2x_2 = 2\) are the same line drawn twice. Every point on it is a solution. Elimination produces a row of all zeros, including on the right: \(0 = 0\), true but uninformative.
The unknowns not attached to a pivot are free variables, and each one adds a dimension to the solution set.
WarningWatch out
R’s solve() raises an error on an exactly singular matrix, but floating-point arithmetic rarely produces an exact zero. A matrix that is singular in theory usually looks nearly singular in practice, and solve() will happily return a confident, meaningless answer. Section 6.10 is about detecting that.
Note what is not on the list: a square system never has exactly two solutions, or seventeen. If two distinct solutions exist, every point on the line between them is also a solution, so there are infinitely many. Linearity does not permit anything in between.
6.6 LU decomposition
Elimination does a lot of work, and all of it depends only on \(\mathbf{A}\) — not on \(\mathbf{b}\). If you have to solve \(\mathbf{A}\mathbf{x} = \mathbf{b}\) for many different \(\mathbf{b}\), you would rather not repeat it.
LU decomposition saves the work. It factors
\[ \mathbf{A} = \mathbf{L}\mathbf{U} \]
where \(\mathbf{L}\) is lower triangular with ones on the diagonal and \(\mathbf{U}\) is upper triangular. \(\mathbf{U}\) is exactly what elimination produced; \(\mathbf{L}\) records the multipliers used along the way.
From the worked example: we subtracted \(2\times\) row 1 from row 2, \(1\times\) row 1 from row 3, then \(2\times\) row 2 from row 3. Those three multipliers are \(\mathbf{L}\):
\[ \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & 2 & 1 \end{bmatrix} \qquad \mathbf{U} = \begin{bmatrix} 1 & 1 & 1 \\ 0 & -1 & -3 \\ 0 & 0 & 7 \end{bmatrix} \]
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 -1
[3,] 1 -1 2[1] TRUEWith the factorization in hand, \(\mathbf{A}\mathbf{x} = \mathbf{b}\) becomes \(\mathbf{L}\mathbf{U}\mathbf{x} = \mathbf{b}\), solved as two triangular systems:
- Forward-substitute \(\mathbf{L}\mathbf{y} = \mathbf{b}\) for \(\mathbf{y}\).
- Back-substitute \(\mathbf{U}\mathbf{x} = \mathbf{y}\) for \(\mathbf{x}\).
Each triangular solve costs \(O(n^2)\), against \(O(n^3)\) for the factorization. Factor once, then every new right-hand side is cheap.
y <- forwardsolve(L, b3)
backsolve(U, y)[1] 1 2 3In practice the factorization includes row swaps for numerical stability, giving \(\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}\) with a permutation matrix \(\mathbf{P}\). That is what LAPACK computes, and what sits behind solve().
NoteIn machine learning
Any time the same matrix is reused against many right-hand sides, factor it once. Gaussian process regression solves against the same kernel matrix repeatedly; iteratively reweighted least squares, behind every glm() fit, refactors a matrix that changes only slightly each step. Recognizing the pattern is worth an order of magnitude.
6.7 Solving systems in R
solve(A3, b3) # solve the system: do this[1] 1 2 3That is the whole API for the common case. What matters is what not to write:
Same answer, more work, worse numerics. Forming \(\mathbf{A}^{-1}\) costs roughly three times a factorization and amplifies rounding error, and you then throw the inverse away. The rule from Section 5.12, stated once more because it is the single most common numerical mistake in applied work: solve, never invert.
Multiple right-hand sides go in as a matrix, and R factors once:
6.8 Overdetermined systems and least squares
Now the case that actually shows up in data science. You have more equations than unknowns — \(n\) observations, \(p\) parameters, \(n \gg p\) — so \(\mathbf{A}\) is tall and thin and the system is overdetermined.
Such a system almost never has a solution. Five points do not lie exactly on a line, and demanding that they do is asking for something false.
So we change the question. Instead of demanding \(\mathbf{X}\boldsymbol{\beta} = \mathbf{y}\) exactly, find the \(\boldsymbol{\beta}\) that comes closest:
\[ \hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 \tag{6.1}\]
Minimizing the squared \(L_2\) norm of the residual is least squares. Squared distance is chosen partly because it is what you want when errors are Gaussian, and partly because it makes the problem solvable in closed form — a fact we lean on throughout.
draw_scatter(x, y, fit = "glm", xlab = "x", ylab = "y")6.9 The normal equations
Equation 6.1 has a closed-form solution, and the geometry of Section 4.14 tells you what it is before any calculus does.
The residual \(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\) is smallest when it is orthogonal to every column of \(\mathbf{X}\) — when there is no part of it left that any predictor could still explain. Orthogonality to every column means
\[ \mathbf{X}^\top(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}) = \mathbf{0} \]
Rearranging gives the normal equations:
\[ \mathbf{X}^\top\mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y} \tag{6.2}\]
This is now a square \(p \times p\) system, and we already know how to solve those.
Worked example. With our five points:
\[ \mathbf{X}^\top\mathbf{X} = \begin{bmatrix} 5 & 15 \\ 15 & 55 \end{bmatrix} \qquad \mathbf{X}^\top\mathbf{y} = \begin{bmatrix} 30.1 \\ 110.2 \end{bmatrix} \]
So the fitted line is \(\hat{y} \approx 0.05 + 1.99x\). The residuals are orthogonal to both columns of \(\mathbf{X}\), as promised:
Two zeros, up to rounding. That is not a coincidence — it is the definition of the solution.
WarningWatch out
Do not compute \(\hat{\boldsymbol{\beta}}\) as solve(t(X) %*% X) %*% t(X) %*% y, even though every textbook writes \((\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\). Forming \(\mathbf{X}^\top\mathbf{X}\) squares the condition number, so a design matrix that was merely awkward becomes numerically hopeless.
Production code uses a QR decomposition of \(\mathbf{X}\) instead, which never forms \(\mathbf{X}^\top\mathbf{X}\) at all. That is what lm() does, and what Section 9.3 covers.
6.10 Conditioning and numerical stability
Here is the failure mode nobody warns you about: a system can have a unique solution that is worthless.
Take two lines that are nearly, but not exactly, parallel.
[1] 0.05[1] 1 1The solution is \((1, 1)\). Now nudge the right-hand side by \(0.05\) — a change of about 2% in one number:
The answer moves from \((1, 1)\) to \((0, 2)\). A tiny change in the data produced a total change in the answer.
The figure explains it. When two lines cross at a shallow angle, the crossing point is poorly pinned down — slide one line a hair and the intersection slides a long way.
The condition number \(\kappa(\mathbf{A})\) measures this. Roughly, it is the factor by which a relative error in \(\mathbf{b}\) can be amplified in \(\mathbf{x}\):
\[ \frac{\text{relative error in } \mathbf{x}}{\text{relative error in } \mathbf{b}} \leq \kappa(\mathbf{A}) \]
Rules of thumb:
| \(\kappa\) | Interpretation |
|---|---|
| near 1 | well-conditioned; errors are not amplified |
| \(10^k\) | expect to lose about \(k\) digits of accuracy |
| \(> 10^{15}\) | nothing survives double precision; treat as singular |
Use \(\kappa\), not \(\det\), to judge a matrix. The determinant has the wrong units: scaling a matrix by \(0.1\) divides the determinant of a \(10 \times 10\) matrix by \(10^{10}\) without making the system any harder. The condition number is scale-invariant.
6.11 Why regularization helps
Ill-conditioning is not a rare numerical curiosity. It is what collinear predictors look like, and collinear predictors are everywhere: height and weight, income and education, any two sensors measuring almost the same thing.
Here are two predictors that are nearly identical, and a response that depends on their sum:
[1] 1114797A condition number around \(10^6\). Least squares still returns an answer:
Look at what it did. It assigned roughly \(-13.6\) to one predictor and \(+15.6\) to the other — enormous, opposite coefficients on two variables that are nearly the same thing. Their sum is about \(1.99\), close to the true value of 2, so the predictions are fine. The individual coefficients are nonsense.
Ridge regression fixes this by adding a small multiple of the identity before solving:
\[ \hat{\boldsymbol{\beta}}_{\text{ridge}} = (\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^\top\mathbf{y} \tag{6.3}\]
Adding \(\lambda\) to the diagonal pulls every eigenvalue up by \(\lambda\) (Section 8.1), which drags the condition number down and makes the problem well-posed again.
ridge <- function(X, y, lambda) {
p <- ncol(X)
pen <- diag(lambda, p)
pen[1, 1] <- 0 # never penalize the intercept
drop(solve(crossprod(X) + pen, crossprod(X, y)))
}
lambdas <- c(0, 0.001, 0.01, 0.1, 1, 10)
paths <- sapply(lambdas, function(l) ridge(Xc, yc, l))
rownames(paths) <- c("intercept", "x1", "x2")
round(paths, 3) [,1] [,2] [,3] [,4] [,5] [,6]
intercept 3.086 3.097 3.104 3.111 3.170 3.729
x1 -13.579 -4.996 0.043 0.892 0.977 0.936
x2 15.569 6.984 1.944 1.093 0.997 0.937By \(\lambda = 1\) the two coefficients are both about \(0.98\) — the model has decided, sensibly, that two indistinguishable predictors should share the effect equally. Their sum stays near 2 throughout, because that part of the problem was always well determined.
NoteIn machine learning
This is the numerical half of why regularization works. The statistical half is the bias–variance tradeoff (Section 21.10): ridge accepts a little bias in exchange for a large drop in variance. Both halves describe the same act — refusing to take the data quite as literally as it asks to be taken.
6.12 Summary
| Task | R |
|---|---|
| Solve a square system | solve(A, b) |
| Multiple right-hand sides | solve(A, B) |
| Triangular solves |
forwardsolve(L, b), backsolve(U, y)
|
| Least squares |
lm(y ~ x) or qr.solve(X, y)
|
| Normal equations (illustration only) | solve(crossprod(X), crossprod(X, y)) |
| Condition number | kappa(A, exact = TRUE) |
| Determinant | det(A) |
6.13 Exercises
1. Solve by hand using elimination, then check in R:
\[ \begin{aligned} 3x_1 + 2x_2 &= 7 \\ x_1 - x_2 &= -1 \end{aligned} \]
TipSolution
2. For each system, say whether it has one solution, none, or infinitely many — without computing anything.
\[ \text{(a)}\;\begin{aligned} x_1 + x_2 &= 3 \\ 2x_1 + 2x_2 &= 6 \end{aligned} \qquad \text{(b)}\;\begin{aligned} x_1 + x_2 &= 3 \\ x_1 + x_2 &= 4 \end{aligned} \qquad \text{(c)}\;\begin{aligned} x_1 + x_2 &= 3 \\ x_1 - x_2 &= 1 \end{aligned} \]
TipSolution
Infinitely many. The second equation is twice the first — one line, drawn twice.
None. Same left-hand side, different right-hand side: parallel lines. The system says \(3 = 4\).
Exactly one. The lines have different slopes, so they cross once, at \((2, 1)\).
3. Verify that \(\mathbf{L}\mathbf{U} = \mathbf{A}_3\) and use the factorization to solve \(\mathbf{A}_3\mathbf{x} = (3, 0, 0)\) without calling solve().
TipSolution
b_new <- c(3, 0, 0)
y_new <- forwardsolve(L, b_new)
x_new <- backsolve(U, y_new)
x_new[1] -0.4285714 2.1428571 1.2857143A3 %*% x_new # should recover b_new [,1]
[1,] 3.000000e+00
[2,] 2.220446e-16
[3,] 1.332268e-15The factorization was computed once and is reused for free on every new right-hand side.
4. Fit a line to x <- 1:5, y <- c(2.1, 3.9, 6.2, 7.8, 10.1) three ways: the normal equations, qr.solve(), and lm(). Confirm they agree.
TipSolution
x
normal 0.05 1.99
qr 0.05 1.99
lm 0.05 1.99They agree to printed precision here because the problem is well-conditioned. On a collinear design they would not, and only the first would be wrong.
5. The Hilbert matrix \(h_{ij} = 1/(i + j - 1)\) is the classic ill-conditioned example. Build the \(6 \times 6\) Hilbert matrix, compute its condition number, then solve \(\mathbf{H}\mathbf{x} = \mathbf{b}\) where \(\mathbf{b}\) is chosen so the exact answer is a vector of ones. How many digits survive?
TipSolution
[1] 14951059[1] 1 1 1 1 1 1The condition number is about \(1.5 \times 10^7\). Double precision carries about 16 digits, so the rule of thumb says up to 7 of them are at risk. The recovered answer is correct to roughly 10 decimal places — comfortably inside that bound, which is the point: \(\kappa\) gives a worst case, not a prediction.
Push further and the bound bites. At \(n = 13\) the condition number exceeds \(10^{18}\), past what double precision can represent at all, and the matrix is numerically singular:
6. Using the collinear design Xc from Section 6.11, compare the predictions from ordinary least squares and from ridge with \(\lambda = 1\). Then compare the coefficients. What do you conclude?
TipSolution
[1] 0.3452561rbind(ols = b_ols, ridge = b_rdg) # coefficients x1 x2
ols 3.085908 -13.5792354 15.5690704
ridge 3.170482 0.9771824 0.9972754The predictions barely differ; the coefficients differ enormously. The data determines what the sum of the two effects must be, but says almost nothing about how to divide it between two nearly identical predictors.
The lesson generalizes: on collinear data, a model can predict well while its coefficients are meaningless. If you intend to interpret coefficients, check the conditioning first.