Solving Systems of Linear Equations

Now that we have matrices, we can formalize how to *solve* $A\mathbf{x} = \mathbf{b}$. The canonical algorithm is Gaussian elimination. The idea: apply a sequence of elementary row operations to reduce $A$ (augmented with $\mathbf{b}$) to an upper-triangular form called row-echelon form, then back-substitute. The three legal operations are swapping two rows, scaling a row by a non-zero constant, and adding a multiple of one row to another.

A stronger target is reduced row-echelon form (RREF): every pivot is 1, and each pivot column has zeros everywhere else. From RREF the solution (or the parametric family of solutions) can be read off directly. Because row operations correspond to multiplying by invertible elementary matrices, they do not change the solution set — just the ease of reading it.

Variables corresponding to pivot columns are basic variables; the others are free variables. When every column has a pivot, the solution is unique. When some columns lack pivots, each free variable contributes an extra dimension to the solution set. When the augmented column has a pivot but the coefficient column does not (a row like $[0 \, 0 \, \dots \, 0 \,|\, 1]$), the system is inconsistent and has no solution.

The general solution of a consistent system is $\mathbf{x} = \mathbf{x}_p + \mathbf{x}_h$, where $\mathbf{x}_p$ is any *particular* solution to $A\mathbf{x} = \mathbf{b}$ and $\mathbf{x}_h$ is any solution to the *homogeneous* system $A\mathbf{x} = \mathbf{0}$. The homogeneous solutions form a vector space — the null space of $A$ — whose dimension equals the number of free variables. This structural view is central to later topics like rank and linear mappings.

Worked example — Gaussian elimination on a 3×3 system: Solve $x + y + z = 6$, $2x - y + z = 3$, $x + 2y - z = 2$. Write the augmented matrix $\left[\begin{smallmatrix}1 & 1 & 1 & | & 6 \\ 2 & -1 & 1 & | & 3 \\ 1 & 2 & -1 & | & 2\end{smallmatrix}\right]$. R2 ← R2 − 2R1 gives row $(0, -3, -1\,|\,-9)$. R3 ← R3 − R1 gives $(0, 1, -2\,|\,-4)$. R3 ← R3 + R2/3 gives $(0, 0, -7/3\,|\,-7)$, so $z = 3$. Back-substitute: $-3y - 3 = -9 \Rightarrow y = 2$, then $x + 2 + 3 = 6 \Rightarrow x = 1$. Solution $(1, 2, 3)$ ✓.

Worked example — inconsistent system detected: Consider $\left[\begin{smallmatrix}1 & 1 & | & 2 \\ 2 & 2 & | & 5\end{smallmatrix}\right]$. R2 ← R2 − 2R1 gives $(0, 0 \,|\, 1)$ — a row that reads $0 = 1$, which is impossible. The system has no solution. Seeing a row of zeros on the left with a non-zero right-hand side is the algebraic fingerprint of inconsistency.

In practice, modern ML code uses numerically stable variants of Gaussian elimination: LU decomposition (factor $A = LU$ into lower and upper triangular matrices) and QR decomposition (factor $A = QR$ into orthogonal and upper-triangular). The LU approach powers most dense direct solvers; QR is more robust for least-squares problems. Under the hood, numpy.linalg.solve and torch.linalg.solve are doing exactly this.