3. Systems of Linear Equations

A system of linear equations asks for the values of several unknowns that satisfy several linear conditions at once. It is the most common computational problem in applied mathematics, behind circuit simulators, structural analysis, economic models, and the inner loops of machine learning. This chapter describes what solution sets of linear systems look like, encodes systems as matrices, and develops Gaussian elimination, the algorithm that solves them. Matrix algebra, determinants, vector spaces, and rank all grow out of the questions first asked here.

<a id="31-linear-systems-and-solutions"></a>

3.1 Linear Systems and Solutions

This subchapter defines linear equations and systems, sorts systems by whether they have solutions, and shows why exactly three outcomes are possible: one solution, infinitely many, or none.

Linear Equations in n Variables

Why it matters. Linear equations are the simplest equations that relate several unknowns: each unknown is multiplied by a constant and the results are added. They model budgets, mixtures, balances, and flows. Nonlinear problems are routinely approximated by linear ones, as in Newton's method or the linearization of a physical model. Knowing exactly what counts as linear tells you when the tools of this chapter apply.

Definition. A linear equation in the variables $x_1, \dots, x_n$ is an equation that can be written as

$$a_1 x_1 + a_2 x_2 + \cdots + a_n x_n = b,$$

where the coefficients $a_1, \dots, a_n$ and the constant term $b$ are numbers from a field, usually $\mathbb{R}$. In dot-product notation (Section 2.3) this is $\mathbf{a} \cdot \mathbf{x} = b$.

A shopping receipt is a good analogy. If $a_j$ is the price of item $j$ and $x_j$ is how many you bought, then $a_1x_1 + \cdots + a_nx_n$ is the total, and the equation says the total equals $b$. Each unknown contributes in proportion to its own amount, and the contributions add up independently. Linearity means exactly this.

What is and is not linear. The variables may appear only to the first power, multiplied by constants, and never inside other functions or multiplied together.

Equation	Linear?	Reason
$3x_1 - 2x_2 + x_3 = 7$	Yes	Standard form
$\sqrt{2},x - \pi y = 0$	Yes	Coefficients can be any constants
$x_1 = 4 - 5x_2$	Yes	Rearranges to $x_1 + 5x_2 = 4$
$xy = 1$	No	Product of two variables
$x^2 + y = 3$	No	Variable raised to a power
$\sin x + y = 0$	No	Variable inside a function
$2/x + y = 1$	No	Variable in a denominator

Solutions. A solution is a vector $(s_1, \dots, s_n)$ that makes the equation true when $x_i = s_i$ for every $i$. The solution set is the set of all solutions. If at least one coefficient is nonzero, the solution set is a hyperplane with normal vector $\mathbf{a}$ (Section 2.5): a line in $\mathbb{R}^2$ and a plane in $\mathbb{R}^3$. You can choose $n - 1$ of the variables freely and solve for the remaining one, as long as that variable has a nonzero coefficient.

Example. For $2x + 3y - z = 6$, solve for $z$: $z = 2x + 3y - 6$. Every solution has the form $(s,\ t,\ 2s + 3t - 6)$ for real numbers $s, t$. The solution set is a plane, and $(1, 1, -1)$ is one of its points.

  y
  |
  2 *                    2x + 3y = 6 in R^2:
  |   `  .               a line with intercepts (3, 0) and (0, 2)
  |        `  .          and normal vector (2, 3)
  |              `  .
  +--------------------*------- x
  0                    3

The degenerate equation. If every $a_i = 0$, the equation reads $0 = b$. When $b = 0$, every vector is a solution. When $b \ne 0$, no vector is. The equation $0 = b$ with $b \ne 0$ will become the signal that a whole system has no solution.

import sympy as sp

x, y, z = sp.symbols("x y z")

def is_linear(expr, variables):
    """True if expr is a polynomial of total degree <= 1 in the given variables."""
    try:
        return sp.Poly(expr, *variables).total_degree() <= 1
    except sp.PolynomialError:                       # sin(x) and 1/x are not polynomials at all
        return False

candidates = [3*x - 2*y + z - 7, sp.sqrt(2)*x - sp.pi*y, x*y - 1,
              x**2 + y - 3, sp.sin(x) + y, 2/x + y - 1]
for e in candidates:
    print(f"{str(e):<22} {is_linear(e, [x, y, z])}")
# True, True, False, False, False, False

eq = sp.Eq(2*x + 3*y - z, 6)
print(sp.solve(eq, z))                               # [2*x + 3*y - 6]: x and y are free
print(eq.subs({x: 1, y: 1, z: -1}))                  # True: (1, 1, -1) is a solution

Key takeaways

• A linear equation has the form $a_1x_1 + \cdots + a_nx_n = b$, or $\mathbf{a} \cdot \mathbf{x} = b$, with constant coefficients.
• Products of variables, powers, roots, and functions of variables make an equation nonlinear.
• A linear equation with a nonzero coefficient vector has a hyperplane as its solution set, with $n - 1$ degrees of freedom.
• The degenerate equation $0 = b$ has no solutions when $b \ne 0$ and has every vector as a solution when $b = 0$.

🧪 Practice

1. Which are linear in $x, y, z$? (a) $4x - y = z + 2$ (b) $x + yz = 1$ (c) $x\cos(\pi/3) + e^2 y = 0$ (d) $\sqrt{x} + y = 3$.
2. Describe all solutions of $x_1 - 4x_2 + 2x_3 = 5$ by solving for $x_1$ in terms of two parameters, and give three specific solutions.
3. For which values of $k$ is the equation $(k^2 - 1)x + (k - 1)y = k + 1$ degenerate, and what is its solution set in each of those cases?
4. Interview: A model predicts a house price as $w_1 \cdot \text{area} + w_2 \cdot \text{rooms} + w_3$. Is the model linear in the features, in the weights, or both? Which view matters when you fit the weights to data? Hint: during fitting, the data values are known numbers and the weights are the unknowns.

Systems of Linear Equations

Why it matters. Real problems impose several conditions at once. A diet must meet several nutrient targets, and a bridge truss must balance forces in each direction at each joint. Each condition is a linear equation, and a solution must satisfy all of them simultaneously.

Definition. A system of $m$ linear equations in $n$ unknowns (an $m \times n$ linear system) is a collection

$$ \begin{aligned} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &= b_1 \ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &= b_2 \ &\ \ \vdots \ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n &= b_m \end{aligned} $$

The double index follows the convention row, then column: $a_{ij}$ is the coefficient of $x_j$ in equation $i$. A solution of the system is a vector that satisfies every equation. The solution set of the system is therefore the intersection of the solution sets of the individual equations. Think of each equation as a filter, and a solution is whatever passes through all of them.

Equivalent systems. Two systems are equivalent if they have the same solution set. Every method for solving systems replaces the system with an equivalent one that is easier to read. Three moves never change the solution set: reordering equations, multiplying an equation by a nonzero constant, and adding a multiple of one equation to another. Section 3.3 explains why.

Example. Solve $$x + 2y = 5, \qquad 3x - y = 1.$$ Subtract 3 times the first equation from the second: $-7y = -14$, so $y = 2$. Substituting back gives $x = 5 - 4 = 1$. The new system ${x + 2y = 5,\ -7y = -14}$ is equivalent to the original but easier to solve.

Triangular systems and back substitution. A system in which each equation involves one fewer unknown than the one above it can be solved from the bottom up:

   x +  y +  z = 6          3z = 9            ->  z = 3
        y + 2z = 8          y + 2(3) = 8      ->  y = 2
            3z = 9          x + 2 + 3 = 6     ->  x = 1
   (triangular shape)       (solve bottom-up, substituting known values)

This process is called back substitution. Gaussian elimination (Section 3.3) is a method for bringing any system into a triangular-like shape so that back substitution can finish the job.

Shapes of systems. The relation between $m$ and $n$ suggests, but does not determine, how many solutions to expect:

Shape	Name	Typical behavior	Exceptions
$m = n$	Square	Exactly one solution	Dependent or contradictory equations
$m < n$	Underdetermined	Infinitely many solutions	Can still have none
$m > n$	Overdetermined	No solution	Can still be consistent if equations are redundant

from fractions import Fraction as F

# Each equation is (coefficients, right-hand side): x + y + z = 6, y + 2z = 8, 3z = 9
system = [([1, 1, 1], 6),
          ([0, 1, 2], 8),
          ([0, 0, 3], 9)]

def satisfies(system, x):
    """A vector solves the system only if it satisfies every equation."""
    return all(sum(a * xi for a, xi in zip(coeffs, x)) == b for coeffs, b in system)

def back_substitute(system):
    """Solve a triangular system: equation i involves only x_i, ..., x_n, with a nonzero coefficient on x_i."""
    n = len(system)
    x = [F(0)] * n
    for i in range(n - 1, -1, -1):                              # last equation first
        coeffs, b = system[i]
        known = sum(coeffs[j] * x[j] for j in range(i + 1, n))  # contribution of already-solved unknowns
        x[i] = (b - known) / F(coeffs[i])
    return x

x = back_substitute(system)
print([str(v) for v in x])            # ['1', '2', '3']
print(satisfies(system, x))           # True
print(satisfies(system, [2, 1, 3]))   # False: equation 1 holds, but 1 + 2*3 = 7 != 8

Key takeaways

• A linear system is a set of linear equations in the same unknowns, and $a_{ij}$ is the coefficient of $x_j$ in equation $i$.
• A solution must satisfy every equation, so the solution set is the intersection of the individual solution sets.
• Solving means replacing the system with an equivalent, simpler one that has the same solution set.
• Triangular systems are solved by back substitution, from the last equation up.

🧪 Practice

1. Solve by elimination: $2x - y = 7$, $x + 3y = 0$.
2. Solve by back substitution: $x_1 - x_2 + 2x_3 = 3$, $\ 2x_2 - x_3 = 1$, $\ 4x_3 = 8$.
3. Give an example of a $3 \times 2$ system that has exactly one solution, and explain why this does not contradict the "typical behavior" column of the table.
4. Interview: You are given 10,000 noisy sensor readings, each a linear equation in 5 unknown calibration constants. Should you expect an exact solution? What would you compute instead? Hint: think about what measurement noise does to an overdetermined system, and which quantity you would minimize.

Consistent and Inconsistent Systems

Why it matters. Before computing a solution, you need to know whether one exists. Inconsistent systems arise from contradictory requirements, from modeling mistakes, and almost always from overdetermined systems built from real measurements. Recognizing inconsistency quickly prevents wasted work and points to the right repair, such as fixing the model or switching to least squares (Section 9.5).

Definition. A system is consistent if it has at least one solution and inconsistent if it has none.

How inconsistency shows up. Inconsistency is exposed when elimination produces the degenerate equation $0 = c$ with $c \ne 0$:

$$ \begin{aligned} x + y &= 2 \ 2x + 2y &= 5 \end{aligned} \quad\xrightarrow{\ \text{eq}_2 - 2,\text{eq}_1\ }\quad \begin{aligned} x + y &= 2 \ 0 &= 1 \end{aligned} $$

Geometrically, the two equations describe parallel lines: same normal vector $(1, 1)$, different offsets.

Hidden contradictions. A system can be inconsistent even when no two of its equations contradict each other:

$$x + y = 1, \qquad y + z = 1, \qquad x + 2y + z = 3.$$

Any two of these equations have common solutions. But adding the first two equations gives $x + 2y + z = 2$, which contradicts the third. This is the triangular-prism configuration of three planes from Section 2.5: each pair of planes meets in a line, but no point lies on all three.

Consistency depends on the right-hand side. With the coefficients fixed, some right-hand sides give consistent systems and others do not. For

$$x + 2y = b_1, \qquad 2x + 4y = b_2,$$

subtracting twice the first equation from the second gives $0 = b_2 - 2b_1$. The system is consistent exactly when $b_2 = 2b_1$. Conditions like this are called compatibility conditions. Section 3.2 restates them as "$\mathbf{b}$ lies in the span of the columns."

An analogy is a client specification that says "the total weight must be under 2 kg" and "each of the two parts must weigh more than 1.5 kg". Every requirement is reasonable on its own, but together they cannot all be met.

Signal during elimination	Meaning
A row $0 = c$ with $c \ne 0$	Inconsistent: stop, no solution exists
A row $0 = 0$	That equation was redundant; drop it
Neither	Consistent; continue to find the solutions

A system whose constant terms are all zero is always consistent, because $\mathbf{x} = \mathbf{0}$ is a solution. These homogeneous systems are studied in Section 3.4.

import sympy as sp

x, y, z = sp.symbols("x y z")

print(sp.linsolve([x + y - 2, 2*x + 2*y - 5], [x, y]))                    # EmptySet
print(sp.linsolve([x + y - 1, y + z - 1, x + 2*y + z - 3], [x, y, z]))    # EmptySet (hidden contradiction)
print(sp.linsolve([x + y - 1, y + z - 1, x + 2*y + z - 2], [x, y, z]))    # {(z, 1 - z, z)}: consistent

# Compatibility condition for x + 2y = b1, 2x + 4y = b2
b1, b2 = sp.symbols("b1 b2")
leftover = sp.expand((2*x + 4*y - b2) - 2*(x + 2*y - b1))   # eq2 - 2*eq1 cancels x and y
print(leftover)                                             # 2*b1 - b2: must be 0 for consistency

Key takeaways

• A system is consistent if it has at least one solution and inconsistent if it has none.
• Elimination exposes inconsistency as an equation $0 = c$ with $c \ne 0$.
• Inconsistency can involve three or more equations even when every pair is compatible.
• For fixed coefficients, consistency depends on the right-hand side through compatibility conditions; homogeneous systems are always consistent.

🧪 Practice

1. Determine whether the system $x - y + z = 2$, $2x + y - z = 1$, $3x + 0y + 0z = 4$ is consistent.
2. For which values of $h$ is the system $x + hy = 4$, $3x + 6y = 8$ consistent? Treat the cases separately.
3. Find the compatibility condition on $b_1, b_2, b_3$ for the system $x + y = b_1$, $y + z = b_2$, $x + 2y + z = b_3$.
4. Interview: A colleague's solver reports "no solution" for a system built from a physical model that must have a solution. List possible causes and how you would diagnose each. Hint: consider modeling errors, units, rounding in near-redundant equations, and whether the solver was even meant for that shape of system.

Geometric Interpretation of Solution Sets

Why it matters. A picture can predict how many solutions to expect before you compute anything. The intuition from lines and planes in two and three dimensions carries over to $\mathbb{R}^n$, where pictures are no longer possible.

The row picture. Each equation $\mathbf{a}_i \cdot \mathbf{x} = b_i$ is a hyperplane: a line in $\mathbb{R}^2$ and a plane in $\mathbb{R}^3$. The solution set of the system is the intersection of these hyperplanes. Two lines in the plane can meet in three ways:

   One point                 No point                  Infinitely many
                             (parallel lines)          (the same line twice)

      \     /                 \      \                      \
       \   /                   \      \                      \\
        \ /                     \      \                      \\
         *                       \      \                      \\
        / \                       \      \                      \
       /   \                    same normal,              equations are
                                different offset          multiples of each other

In $\mathbb{R}^3$, three planes can meet in a point, along a line, in a whole plane, or not at all, with the configurations listed in Section 2.5.

The column picture. The same system can be read differently. Instead of intersecting rows, ask which combination of the coefficient columns produces the right-hand side. Section 3.2 develops this view. The two pictures always give the same answer, and each one makes some questions easier.

Possible shapes. In $\mathbb{R}^n$, the solution set of a linear system is always one of the following:

Solution set	Dimension	Example in $\mathbb{R}^3$
Empty	-	Two parallel planes
A single point	0	Three planes meeting at a corner
A line	1	Two non-parallel planes
A plane	2	One equation, or several copies of the same plane
A $k$-dimensional flat	$k$	$\mathbf{p} + t_1\mathbf{v}_1 + \cdots + t_k\mathbf{v}_k$

A solution set is never curved and never a finite set with two or more points. It is always flat. The next topic explains why.

Counting dimensions. As a rule of thumb, each independent equation removes one dimension. With $n$ unknowns and $m$ independent, consistent equations, the solution set has dimension $n - m$. Two planes in $\mathbb{R}^3$ give $3 - 2 = 1$, a line. The word "independent" carries the weight here, and rank (Section 3.4) makes it precise.

import sympy as sp

x, y, z = sp.symbols("x y z")

def describe(eqs, unknowns):
    """Name the shape of the solution set by counting free parameters."""
    sol = sp.linsolve(eqs, unknowns)
    if sol == sp.EmptySet:
        return "empty"
    (point,) = sol                                         # one parametric tuple
    free = set().union(*(c.free_symbols for c in point)) & set(unknowns)
    shape = {0: "point", 1: "line", 2: "plane", 3: "all of R^3"}[len(free)]
    return f"{shape}: {point}"

print(describe([x + y + z - 6, x - y + z - 2, x - z], [x, y, z]))    # point: (2, 2, 2)
print(describe([x + y + z - 6, x - y + z - 2], [x, y, z]))           # line: (4 - z, 2, z)
print(describe([x + y + z - 6, 2*x + 2*y + 2*z - 12], [x, y, z]))    # plane: (-y - z + 6, y, z)
print(describe([x + y + z - 6, x + y + z - 7], [x, y, z]))           # empty

Key takeaways

• Row picture: each equation is a hyperplane, and the solution set is their intersection.
• Column picture: solving means finding weights that combine the columns into the right-hand side.
• Solution sets are always flat: empty, a point, a line, a plane, or a higher-dimensional flat.
• Each independent equation typically removes one dimension, so $m$ independent equations in $n$ unknowns leave $n - m$ dimensions.

🧪 Practice

1. Sketch the lines $x - y = 0$, $x + y = 2$, and $x = 1$ in the plane. Do all three meet at a common point?
2. Describe geometrically the solution set of the single equation $x + y + z + w = 1$ in $\mathbb{R}^4$. What is its dimension?
3. Give three equations in $x, y, z$ whose solution set is exactly the line through the origin with direction $(1, 2, 3)$. What is the smallest number of equations that can do this?
4. Interview: Explain why two linear equations in three unknowns can never have exactly one solution, using only geometry. Hint: intersect two planes in $\mathbb{R}^3$ and list every possible outcome.

Unique, Infinite, and No Solutions

Why it matters. A linear system has either no solutions, exactly one, or infinitely many. It can never have exactly two, unlike $x^2 = 1$. This trichotomy is one of the most useful facts in the subject. As soon as you find two different solutions, you know there are infinitely many, and the uniqueness question reduces to a question that does not involve the right-hand side at all.

Theorem. Every system of linear equations over $\mathbb{R}$ (or $\mathbb{C}$) has either no solution, exactly one solution, or infinitely many solutions.

Proof. Assume there are two distinct solutions $\mathbf{p} \ne \mathbf{q}$. For any real $t$, consider $\mathbf{r}(t) = (1 - t)\mathbf{p} + t\mathbf{q}$, a point on the line through $\mathbf{p}$ and $\mathbf{q}$. For every equation $\mathbf{a}_i \cdot \mathbf{x} = b_i$ of the system:

$$\mathbf{a}_i \cdot \mathbf{r}(t) = (1 - t),\mathbf{a}_i \cdot \mathbf{p} + t,\mathbf{a}_i \cdot \mathbf{q} = (1 - t)b_i + t b_i = b_i.$$

So every $\mathbf{r}(t)$ is a solution. Because $\mathbf{q} - \mathbf{p} \ne \mathbf{0}$, different values of $t$ give different points, so there are infinitely many solutions. $\blacksquare$

The proof shows more than the count: the entire line through any two solutions consists of solutions, which is why solution sets are flat.

Contrast with nonlinear systems.

System	Number of solutions
$x^2 = 1$	2: $x = \pm 1$
$x^3 - x = 0$	3: $x = -1, 0, 1$
$xy = 0,\ x + y = 1$	2: $(1, 0)$ and $(0, 1)$
Any linear system	0, 1, or infinitely many

Over finite fields. The proof uses infinitely many values of $t$. Over a finite field $\mathbb{F}_q$ (Section 1.2), there are only $q$ choices of $t$. There, the number of solutions is $0$ or $q^k$ for some $k \ge 0$. A system over $\mathbb{F}_2$ can have exactly two solutions ($k = 1$), but never exactly three.

Telling the cases apart. Section 3.3 reduces the system to echelon form, and then:

Echelon form contains	Outcome
A row $[,0 \ \cdots \ 0 \mid c,]$ with $c \ne 0$	No solution
No such row, and every variable is basic	Exactly one solution
No such row, and at least one free variable	Infinitely many solutions

Unique in theory, fragile in practice. Two nearly parallel lines have a unique intersection point, but that point moves a long way when the data change slightly:

$$x + y = 2,\quad x + 1.0001y = 2.0001 \ \Rightarrow\ (1, 1); \qquad x + y = 2,\quad x + 1.0001y = 2.0002 \ \Rightarrow\ (0, 2).$$

A change of $0.0001$ in one constant moves the solution by a distance of about $1.4$. Such systems are called ill-conditioned (Section 11.3).

import numpy as np

# x + 2y - z = 3 and x - y + z = 0, two known solutions p and q
eqs = [(np.array([1., 2., -1.]), 3.0), (np.array([1., -1., 1.]), 0.0)]
p = np.array([1., 1., 0.])
q = np.array([0., 3., 3.])
for t in [-2.0, 0.5, 10.0]:
    r = (1 - t) * p + t * q                               # a point on the line through p and q
    print(t, [float(a @ r) for a, _ in eqs])              # always [3.0, 0.0]: r is a solution too

# Nearly parallel lines: a tiny change in b moves the solution far
A = np.array([[1., 1.], [1., 1.0001]])
print(np.linalg.solve(A, [2., 2.0001]))                   # approximately [1. 1.]
print(np.linalg.solve(A, [2., 2.0002]))                   # approximately [0. 2.]

Key takeaways

• Over $\mathbb{R}$ or $\mathbb{C}$, a linear system has 0, 1, or infinitely many solutions, never exactly 2.
• If $\mathbf{p}$ and $\mathbf{q}$ are solutions, so is every point on the line through them.
• Over a finite field $\mathbb{F}_q$, the count is 0 or a power of $q$.
• A unique solution can still be numerically fragile when the hyperplanes are nearly parallel.

🧪 Practice

1. A linear system is known to have the solutions $(1, 2, 0)$ and $(3, 0, 2)$. Write down three more solutions.
2. Show by example that the system $x^2 + y^2 = 1$, $x = y$ has exactly two solutions. Which step of the proof above fails for it?
3. Over $\mathbb{F}_2$, list all solutions of $x + y + z = 1$ and verify that their number is a power of 2.
4. Interview: A numerical solver returns a unique solution for a $2 \times 2$ system, but changing the fifth decimal of one input changes the answer completely. What is going on, and how would you quantify it? Hint: look at the angle between the two lines, and at the ratio of the change in the output to the change in the input.

<a id="32-matrix-representation"></a>

3.2 Matrix Representation

Writing out variable names in every equation is redundant, so this subchapter reduces a system to its numbers and shows three equivalent ways to read them: as an augmented matrix, as a matrix-vector equation $A\mathbf{x} = \mathbf{b}$, and as a vector equation about linear combinations.

Coefficient Matrix

Why it matters. The names of the unknowns carry no information. Only the coefficients and their positions do. Dropping the names leaves a compact rectangular table that algorithms and computers can process directly. The table also exposes structure, such as many zeros or a symmetric pattern, that the written equations hide.

Definition. The coefficient matrix of an $m \times n$ system is the $m \times n$ array $A = [a_{ij}]$ whose entry in row $i$ and column $j$ is the coefficient of $x_j$ in equation $i$:

• Rows correspond to equations. Row $i$ is the coefficient vector $\mathbf{a}_i$, the normal vector of the $i$-th hyperplane.
• Columns correspond to unknowns. Column $j$ records how $x_j$ enters each equation.
• A variable missing from an equation contributes a 0 entry.

The size is written $m \times n$, rows first. This chapter uses a running example:

$$ \begin{aligned} 2x + y - z &= 8 \ -3x - y + 2z &= -11 \ -2x + y + 2z &= -3 \end{aligned} \qquad\Longrightarrow\qquad A = \begin{bmatrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{bmatrix}. $$

               x    y    z         <- one column per unknown, in a fixed order
            [  2    1   -1 ]       <- equation 1
     A  =   [ -3   -1    2 ]       <- equation 2
            [ -2    1    2 ]       <- equation 3
                    ^
                    a_32 = 1: coefficient of y (column 2) in equation 3

Writing it correctly. Most mistakes come from copying, not from algebra:

Pitfall	Example	Fix
Forgetting a missing variable	$x + 3z = 4$ written as row $[1\ \ 3]$	Insert a zero: $[1\ \ 0\ \ 3]$
Inconsistent variable order	Row 1 in order $x, y, z$; row 2 in order $y, x, z$	Fix one order for all rows
Terms on the wrong side	$2x = 3y + 1$	Move all variables left: $2x - 3y = 1$

Why the shape matters. Large systems from physics or networks often have coefficient matrices that are almost entirely zeros (sparse). Each equation involves only a few neighboring unknowns. Specialized algorithms (Section 12.2) exploit this, and seeing the matrix is how you notice it.

import sympy as sp

x, y, z = sp.symbols("x y z")
eqs = [sp.Eq(2*x + y - z, 8),
       sp.Eq(-3*x - y + 2*z, -11),
       sp.Eq(-2*x + y + 2*z, -3)]

A, b = sp.linear_eq_to_matrix(eqs, [x, y, z])     # the variable order fixes the column order
print(A)          # Matrix([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]])
print(A.shape)    # (3, 3)
print(A[2, 1])    # 1: coefficient of y in equation 3 (Python indices start at 0)

# Missing variables become zeros; terms on the right-hand side are moved automatically
A2, b2 = sp.linear_eq_to_matrix([sp.Eq(x + 3*z, 4), sp.Eq(y, z + 1)], [x, y, z])
print(A2, b2)     # Matrix([[1, 0, 3], [0, 1, -1]]) Matrix([[4], [1]])

Key takeaways

• The coefficient matrix $A = [a_{ij}]$ has one row per equation and one column per unknown.
• Entry $a_{ij}$ is the coefficient of $x_j$ in equation $i$; a missing variable gives a 0.
• All rows must use the same variable order, with every variable term moved to the left side.
• Row $i$ is the normal vector of equation $i$; column $j$ shows how $x_j$ affects every equation.

🧪 Practice

1. Write the coefficient matrix of $x_1 - 2x_3 = 4$, $\ 3x_2 + x_1 = 0$, $\ x_3 = 5 - x_2$, using the order $x_1, x_2, x_3$.
2. Write a system of equations whose coefficient matrix is $\begin{bmatrix} 0 & 1 & -1 & 2 \ 4 & 0 & 0 & 1 \end{bmatrix}$. How many equations and unknowns does it have?
3. If you reorder the unknowns as $z, x, y$, how does the coefficient matrix of the running example change? Does the solution change?
4. Interview: A finite-difference discretization of a 1D heat equation on 1,000 grid points produces a $1000 \times 1000$ coefficient matrix in which each row has at most 3 nonzero entries. Why does that matter for storage and for solving? Hint: compare the number of stored entries in dense and sparse storage, and think about which entries elimination actually touches.

Augmented Matrix

Why it matters. The coefficient matrix alone loses the right-hand side, but elimination must change both sides of each equation together. The augmented matrix keeps them together, and Gaussian elimination operates on it.

Definition. The augmented matrix of the system with coefficient matrix $A$ and right-hand side $\mathbf{b}$ is the $m \times (n+1)$ matrix

$$[,A \mid \mathbf{b},] = \left[\begin{array}{ccc|c} a_{11} & \cdots & a_{1n} & b_1 \ \vdots & & \vdots & \vdots \ a_{m1} & \cdots & a_{mn} & b_m \end{array}\right].$$

The vertical bar marks where the equals sign was. It is only a visual aid; mathematically this is an ordinary matrix with one extra column.

    2x +  y -  z =   8            [  2   1  -1 |   8 ]
   -3x -  y + 2z = -11    <-->    [ -3  -1   2 | -11 ]
   -2x +  y + 2z =  -3            [ -2   1   2 |  -3 ]

   each row is one complete equation, so an operation on a row
   is an operation on that equation (both sides at once)

Reading an augmented matrix back. Each row is an equation. Some rows have special meaning:

Row	Equation	Meaning
$[,1\ \ 0\ \ 0 \mid 2,]$	$x = 2$	A variable is determined
$[,0\ \ 1\ \ 3 \mid 4,]$	$y + 3z = 4$	A relation between variables
$[,0\ \ 0\ \ 0 \mid 0,]$	$0 = 0$	Redundant; carries no information
$[,0\ \ 0\ \ 0 \mid 5,]$	$0 = 5$	Contradiction: the system is inconsistent

Several right-hand sides at once. If you need to solve $A\mathbf{x} = \mathbf{b}_1$ and $A\mathbf{x} = \mathbf{b}_2$ with the same $A$, augment with both columns: $[,A \mid \mathbf{b}_1\ \mathbf{b}_2,]$. Elimination depends only on $A$, so one pass handles both. Section 4.3 uses this to compute inverses with $[,A \mid I,]$.

import numpy as np

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
b = np.array([8., -11, -3])
M = np.column_stack([A, b])                # augmented matrix [A | b], shape (3, 4)
print(M.shape)                             # (3, 4)

def row_to_equation(row, names="xyz"):
    terms = " ".join(f"{c:+g}{v}" for c, v in zip(row[:-1], names) if c != 0)
    return f"{terms or '0'} = {row[-1]:g}"  # an all-zero coefficient row reads "0 = c"

for row in M:
    print(row_to_equation(row))            # +2x +1y -1z = 8, and so on
print(row_to_equation(np.array([0., 0, 0, 5])))   # 0 = 5: the sign of an inconsistent system

# Two right-hand sides share one coefficient matrix: [A | b1 b2]
b2 = np.array([1., 0, 0])
print(np.column_stack([A, b, b2]).shape)   # (3, 5)

Key takeaways

• The augmented matrix $[,A \mid \mathbf{b},]$ stores the entire system, coefficients and constants, in $m \times (n + 1)$ numbers.
• Each row is one equation, so row operations are operations on equations.
• A zero row $[,0 \cdots 0 \mid 0,]$ is redundant; $[,0 \cdots 0 \mid c,]$ with $c \ne 0$ means no solution.
• Augmenting with several right-hand sides solves several systems with the same $A$ in one elimination.

🧪 Practice

1. Write the augmented matrix of $x - y = 3$, $\ 2y + z = -1$, $\ x + z = 0$.
2. Write the system represented by $\left[\begin{array}{ccc|c} 1 & 0 & -2 & 5 \ 0 & 1 & 4 & -1 \ 0 & 0 & 0 & 0 \end{array}\right]$, and describe its solutions.
3. An augmented matrix has 4 rows and 6 columns. How many equations and unknowns does the system have?
4. Interview: You must solve $A\mathbf{x} = \mathbf{b}$ for the same $1000 \times 1000$ matrix $A$ and 500 different vectors $\mathbf{b}$ that arrive one at a time. Why is running elimination 500 times wasteful, and what would you do instead? Hint: the row operations depend only on $A$; think about recording them once (Section 4.4).

Matrix-Vector Form Ax = b

Why it matters. The form $A\mathbf{x} = \mathbf{b}$ compresses a system of any size into one equation that looks like the scalar equation $ax = b$. The resemblance suggests ideas such as "dividing by $A$", which leads to the inverse in Section 4.3. It also lets software solve a system with one call. It presents $A$ as a function $\mathbf{x} \mapsto A\mathbf{x}$, the starting point of linear transformations (Chapter 7).

The matrix-vector product. For an $m \times n$ matrix $A$ and a vector $\mathbf{x} \in \mathbb{R}^n$, the product $A\mathbf{x}$ is the vector in $\mathbb{R}^m$ whose $i$-th entry is the dot product of row $i$ of $A$ with $\mathbf{x}$:

$$(A\mathbf{x})i = a{i1}x_1 + a_{i2}x_2 + \cdots + a_{in}x_n = \mathbf{a}_i \cdot \mathbf{x}.$$

The product is defined only when the number of columns of $A$ equals the number of entries of $\mathbf{x}$:

        A          x          b
     (m x n)    (n x 1)    (m x 1)
          \______/
       these must match;  the result has m entries, one per equation

With this product, the whole system reads $A\mathbf{x} = \mathbf{b}$: "row $i$ of $A$ dotted with $\mathbf{x}$ equals $b_i$, for every $i$" is exactly the $i$-th equation.

Example. For the running example and $\mathbf{x} = (2, 3, -1)$:

$$ \begin{bmatrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{bmatrix} \begin{bmatrix} 2 \ 3 \ -1 \end{bmatrix}

\begin{bmatrix} 4 + 3 + 1 \ -6 - 3 - 2 \ -4 + 3 - 2 \end{bmatrix}

\begin{bmatrix} 8 \ -11 \ -3 \end{bmatrix} = \mathbf{b}, $$

so $(2, 3, -1)$ solves the system.

Linearity. For all vectors $\mathbf{u}, \mathbf{v}$ and scalars $c$:

$$A(\mathbf{u} + \mathbf{v}) = A\mathbf{u} + A\mathbf{v}, \qquad A(c\mathbf{u}) = c,(A\mathbf{u}).$$

Both follow from the corresponding properties of the dot product. These two rules explain the flatness of solution sets (Section 3.1) and the structure theorems of Section 3.4.

The function view. Think of $A$ as a machine that takes an input $\mathbf{x} \in \mathbb{R}^n$ and outputs $A\mathbf{x} \in \mathbb{R}^m$. If $\mathbf{x}$ lists the amounts of ingredients in a recipe and row $i$ of $A$ lists the amount of nutrient $i$ per unit of each ingredient, then $A\mathbf{x}$ lists the total nutrients. Computing $A\mathbf{x}$ is the forward problem. Solving $A\mathbf{x} = \mathbf{b}$ is the inverse problem: find the ingredient amounts that hit the nutrient targets $\mathbf{b}$.

import numpy as np

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
b = np.array([8., -11, -3])
x = np.array([2., 3, -1])

row_view = np.array([A[i] @ x for i in range(A.shape[0])])  # entry i = (row i of A) . x
print(row_view, A @ x)                                      # both [  8. -11.  -3.]
print(np.allclose(A @ x, b))                                # True: x solves the system

print(np.linalg.solve(A, b))       # approximately [ 2.  3. -1.] (square systems with a unique solution only)

# Linearity: A(u + v) = Au + Av and A(cu) = c(Au)
u, v = np.array([1., 0, 2]), np.array([-1., 4, 0])
print(np.allclose(A @ (u + v), A @ u + A @ v), np.allclose(A @ (3 * u), 3 * (A @ u)))   # True True

Key takeaways

• $A\mathbf{x}$ is the vector whose $i$-th entry is (row $i$ of $A$) $\cdot\ \mathbf{x}$; it requires $A$ to have as many columns as $\mathbf{x}$ has entries.
• The equation $A\mathbf{x} = \mathbf{b}$ is the whole system written as one equation.
• Matrix-vector multiplication is linear: $A(\mathbf{u} + \mathbf{v}) = A\mathbf{u} + A\mathbf{v}$ and $A(c\mathbf{u}) = cA\mathbf{u}$.
• Solving $A\mathbf{x} = \mathbf{b}$ means finding the input that the map $\mathbf{x} \mapsto A\mathbf{x}$ sends to $\mathbf{b}$.

🧪 Practice

1. Compute $\begin{bmatrix} 1 & -2 & 0 \ 3 & 1 & 4 \end{bmatrix}\begin{bmatrix} 2 \ 1 \ -1 \end{bmatrix}$, and state the system of equations that this product turns into when the result is set equal to $(0, 3)$.
2. Explain why $\begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix}\begin{bmatrix} 1 \ 2 \ 3 \end{bmatrix}$ is undefined.
3. Suppose $A\mathbf{u} = (1, 2)$ and $A\mathbf{v} = (0, 5)$. Compute $A(2\mathbf{u} - \mathbf{v})$ without knowing $A$.
4. Interview: In a neural network, a dense layer computes $W\mathbf{x} + \mathbf{c}$. Given the layer's output and $W$, when can you recover the input $\mathbf{x}$ uniquely? Hint: rewrite the problem as a linear system and compare the number of unknowns with the number of independent equations.

Vector Equation Form

Why it matters. The same system can be read column by column. It asks which weights combine the column vectors of $A$ into $\mathbf{b}$. This column picture turns solvability into a geometric question about spans (Section 2.2), and it leads directly to column spaces and linear independence (Chapter 6).

The column view of $A\mathbf{x}$. Write $A$ in terms of its columns, $A = [,\mathbf{c}_1\ \ \mathbf{c}_2\ \ \cdots\ \ \mathbf{c}_n,]$. Then

$$A\mathbf{x} = x_1\mathbf{c}_1 + x_2\mathbf{c}_2 + \cdots + x_n\mathbf{c}_n.$$

The $i$-th entry of the right-hand side is $x_1a_{i1} + \cdots + x_na_{in}$, which matches the row view entry by entry. So the system is equivalent to the vector equation

$$x_1\mathbf{c}_1 + x_2\mathbf{c}_2 + \cdots + x_n\mathbf{c}_n = \mathbf{b}.$$

Example. The system $x - y = 1$, $x + y = 3$ becomes

$$x\begin{bmatrix} 1 \ 1 \end{bmatrix} + y\begin{bmatrix} -1 \ 1 \end{bmatrix} = \begin{bmatrix} 1 \ 3 \end{bmatrix}.$$

The solution $x = 2$, $y = 1$ says: two copies of $\mathbf{c}_1$ followed by one copy of $\mathbf{c}_2$ land exactly on $\mathbf{b}$.

      y
    3 |   * b = (1, 3)
      |     \
    2 |       * 2c1 = (2, 2)         b = 2 c1 + 1 c2
      |     /                        c1 = (1, 1),  c2 = (-1, 1)
    1 |   * c1                       (walk along c1 twice, then along c2 once)
      | /
    0 +---+---+---+-- x
      0   1   2   3

Solvability as a span question. The vector equation has a solution exactly when $\mathbf{b}$ is a linear combination of the columns:

$$A\mathbf{x} = \mathbf{b} \text{ is consistent} \iff \mathbf{b} \in \operatorname{Span}{\mathbf{c}_1, \dots, \mathbf{c}_n}.$$

For example, the columns of $\begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}$ both lie on the line $y = 2x$, so their span is that line. The system is solvable for $\mathbf{b} = (1, 2)$ but not for $\mathbf{b} = (1, 1)$. This is the compatibility condition $b_2 = 2b_1$ from Section 3.1, seen geometrically.

Four views of one system.

View	Notation	Question it asks
System of equations	$a_{i1}x_1 + \cdots + a_{in}x_n = b_i$	Which $\mathbf{x}$ satisfies all equations?
Augmented matrix	$[,A \mid \mathbf{b},]$	(The data structure elimination works on)
Matrix equation	$A\mathbf{x} = \mathbf{b}$	Which input does the map $A$ send to $\mathbf{b}$?
Vector equation	$x_1\mathbf{c}_1 + \cdots + x_n\mathbf{c}_n = \mathbf{b}$	Which weights combine the columns into $\mathbf{b}$?

	Row picture	Column picture
Objects	$m$ hyperplanes in $\mathbb{R}^n$	$n$ column vectors in $\mathbb{R}^m$
A solution is	A point on every hyperplane	A set of weights that reaches $\mathbf{b}$
No solution means	The hyperplanes have no common point	$\mathbf{b}$ is outside the span of the columns

import numpy as np

c1, c2 = np.array([1., 1]), np.array([-1., 1])
b = np.array([1., 3])
A = np.column_stack([c1, c2])                 # the columns of A are c1 and c2
w = np.linalg.solve(A, b)
print(w, w[0] * c1 + w[1] * c2)               # [2. 1.] [1. 3.]: b is a combination of the columns

# Column view of the product: A @ x equals the weighted sum of the columns
A3 = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
x3 = np.array([2., 3, -1])
combo = sum(x3[j] * A3[:, j] for j in range(A3.shape[1]))
print(combo, A3 @ x3)                         # identical: [  8. -11.  -3.]

# Is b in the span of the columns? Compare ranks of A and [A | b] (made precise in Section 3.4)
C = np.array([[1., 2], [2, 4]])
for rhs in (np.array([1., 2]), np.array([1., 1])):
    in_span = np.linalg.matrix_rank(C) == np.linalg.matrix_rank(np.column_stack([C, rhs]))
    print(rhs, in_span)                       # [1. 2.] True, then [1. 1.] False

Key takeaways

• $A\mathbf{x} = x_1\mathbf{c}_1 + \cdots + x_n\mathbf{c}_n$: the product is a combination of the columns, weighted by the entries of $\mathbf{x}$.
• The system $A\mathbf{x} = \mathbf{b}$ is equivalent to the vector equation $x_1\mathbf{c}_1 + \cdots + x_n\mathbf{c}_n = \mathbf{b}$.
• $A\mathbf{x} = \mathbf{b}$ is consistent exactly when $\mathbf{b}$ lies in the span of the columns of $A$.
• The row picture (intersecting hyperplanes) and the column picture (combining columns) always agree; use whichever makes the question easier.

🧪 Practice

1. Write $x_1 + 2x_2 - x_3 = 4$, $\ 3x_1 - x_3 = 1$ as a vector equation.
2. Is $(3, 1, 5)$ a linear combination of $(1, 0, 1)$ and $(1, 1, 2)$? Answer by setting up and solving a linear system.
3. Let $\mathbf{c}_1 = (1, 2, 3)$ and $\mathbf{c}_2 = (2, 4, 6)$. Describe geometrically all $\mathbf{b} \in \mathbb{R}^3$ for which $x_1\mathbf{c}_1 + x_2\mathbf{c}_2 = \mathbf{b}$ is consistent.
4. Interview: A $3 \times 3$ matrix has columns that all lie in one plane through the origin. Explain, using the column picture, why $A\mathbf{x} = \mathbf{b}$ cannot be solvable for every $\mathbf{b} \in \mathbb{R}^3$. Hint: every combination of vectors in a plane through the origin stays in that plane.

<a id="33-gaussian-elimination"></a>

3.3 Gaussian Elimination

Gaussian elimination is a systematic version of the "eliminate a variable" technique from school algebra; this subchapter defines the moves it may use, the staircase shapes it aims for, and the algorithm itself.

Elementary Row Operations

Why it matters. To solve a system, you need moves that simplify it without changing its solution set. Exactly three such moves are enough to solve every linear system. Because they act on rows of the augmented matrix, they are mechanical enough for a computer to perform.

The three operations.

Operation	Notation	Effect on equations	Undone by
Interchange	$R_i \leftrightarrow R_j$	Swap two equations	The same swap
Scaling	$R_i \leftarrow cR_i$, $\ c \ne 0$	Multiply both sides of an equation by $c$	$R_i \leftarrow \tfrac{1}{c}R_i$
Replacement	$R_i \leftarrow R_i + cR_j$, $\ i \ne j$	Add a multiple of another equation	$R_i \leftarrow R_i - cR_j$

Two matrices are row equivalent if a sequence of elementary row operations turns one into the other.

Why the solution set does not change. Consider one operation, taking system $S$ to system $S'$.

1. Every equation of $S'$ is a combination of equations of $S$. So any vector satisfying $S$ also satisfies $S'$.
2. The operation can be undone by another elementary operation, which takes $S'$ back to $S$. By the same argument, any vector satisfying $S'$ satisfies $S$.

Therefore $S$ and $S'$ have the same solution set. Think of a balance scale: adding equal weights to both pans, or doubling everything on both pans, keeps it balanced. Each step can be reversed, so no information is lost.

Why the restrictions. Scaling by $c = 0$ turns an equation into $0 = 0$ and cannot be undone, so the system may gain solutions. A "replacement" $R_i \leftarrow R_i - R_i$ with $i = j$ has the same problem. Both restrictions keep every operation reversible.

Example. On the running example, two replacements clear $x$ from equations 2 and 3:

 [  2   1  -1 |   8 ]   R2 <- R2 + (3/2) R1    [ 2   1     -1   |  8 ]
 [ -3  -1   2 | -11 ]  ------------------->    [ 0   1/2   1/2  |  1 ]
 [ -2   1   2 |  -3 ]   R3 <- R3 + R1          [ 0   2     1    |  5 ]

Each operation on rows is equivalent to multiplying by an elementary matrix on the left. Section 4.4 uses this fact to turn elimination into the LU factorization.

from fractions import Fraction as F

def swap(M, i, j):
    M[i], M[j] = M[j], M[i]

def scale(M, i, c):
    if c == 0:
        raise ValueError("scaling by 0 is not reversible")   # it would erase an equation
    M[i] = [c * a for a in M[i]]

def replace(M, i, j, c):
    """R_i <- R_i + c * R_j."""
    if i == j:
        raise ValueError("replacement needs two different rows")
    M[i] = [a + c * aj for a, aj in zip(M[i], M[j])]

M = [[F(v) for v in row] for row in [[2, 1, -1, 8], [-3, -1, 2, -11], [-2, 1, 2, -3]]]
original = [row[:] for row in M]

replace(M, 1, 0, F(3, 2))     # R2 <- R2 + (3/2) R1 clears x from equation 2
replace(M, 2, 0, F(1))        # R3 <- R3 + R1 clears x from equation 3
for row in M:
    print([str(a) for a in row])
# ['2', '1', '-1', '8']
# ['0', '1/2', '1/2', '1']
# ['0', '2', '1', '5']

replace(M, 2, 0, F(-1))       # undo, in reverse order, with the opposite multiples
replace(M, 1, 0, F(-3, 2))
print(M == original)          # True

Key takeaways

• The elementary row operations are interchange, scaling by a nonzero constant, and replacement $R_i \leftarrow R_i + cR_j$ with $i \ne j$.
• Each operation is reversible by an operation of the same type.
• Reversibility guarantees that row-equivalent augmented matrices describe systems with the same solution set.
• Scaling by zero and adding a row to itself are excluded because they lose information.

🧪 Practice

1. Apply $R_1 \leftrightarrow R_2$, then $R_2 \leftarrow R_2 - 3R_1$ to $\left[\begin{array}{cc|c} 3 & 1 & 2 \ 1 & -1 & 0 \end{array}\right]$. What single operation would undo the second step?
2. Show that $R_1 \leftrightarrow R_2$ can be performed using only replacements and one scaling by $-1$.
3. Explain why the operation $R_1 \leftarrow 2R_1 + R_2$ preserves the solution set even though it is not on the list. Write it as a combination of listed operations.
4. Interview: Is the operation "replace $R_1$ with $R_1 + R_2$ and, simultaneously, $R_2$ with $R_1 + R_2$" safe? Hint: look at the two resulting rows and ask whether the original system can be recovered from them.

Row Echelon Form

Why it matters. Elimination needs a target. Row echelon form is a staircase shape in which back substitution is immediate and consistency can be read off by looking at the last column.

Definition. A matrix is in row echelon form (REF) if

1. all zero rows are at the bottom, and
2. the leading entry of each nonzero row (its first nonzero entry from the left) is strictly to the right of the leading entry of the row above it.

As a consequence, every entry below a leading entry is zero. The leading entries are the pivots of the echelon form.

   [ #  *  *  *  * | * ]       # = leading entry (nonzero)
   [ 0  #  *  *  * | * ]       * = any value
   [ 0  0  0  #  * | * ]
   [ 0  0  0  0  0 | * ]       the steps may be wider than one column,
                               but they always descend to the right

Examples and non-examples.

Matrix	REF?	Reason
$\begin{bmatrix} 2 & 1 & -1 \ 0 & \tfrac12 & \tfrac12 \ 0 & 0 & -1 \end{bmatrix}$	Yes	Leading entries step right in each row
$\begin{bmatrix} 1 & 4 & 0 & 2 \ 0 & 0 & 3 & 1 \end{bmatrix}$	Yes	Steps can skip columns
$\begin{bmatrix} 1 & 2 & 3 \ 0 & 0 & 0 \ 0 & 0 & 4 \end{bmatrix}$	No	Zero row above a nonzero row
$\begin{bmatrix} 1 & 2 & 3 \ 0 & 5 & 1 \ 0 & 3 & 2 \end{bmatrix}$	No	Two leading entries in column 2

Reading an echelon form. For an augmented matrix in REF:

• The system is inconsistent exactly when some row is $[,0 \ \cdots\ 0 \mid c,]$ with $c \ne 0$, that is, when a leading entry sits in the last column.
• Otherwise, back substitution from the bottom row up gives the solutions.

Example. The running example reduces to

$$\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \ 0 & \tfrac12 & \tfrac12 & 1 \ 0 & 0 & -1 & 1 \end{array}\right] \ \Rightarrow\ z = -1,\quad \tfrac12 y - \tfrac12 = 1 \Rightarrow y = 3,\quad 2x + 3 + 1 = 8 \Rightarrow x = 2.$$

Not unique. A matrix has many echelon forms. Scaling a row of an REF gives another REF, for example. The positions of the leading entries, however, are the same in every echelon form of a given matrix (see Pivot Positions below).

def leading_index(row):
    """Column of the first nonzero entry, or None for a zero row."""
    return next((j for j, a in enumerate(row) if a != 0), None)

def is_ref(M):
    last = -1
    seen_zero_row = False
    for row in M:
        lead = leading_index(row)
        if lead is None:
            seen_zero_row = True
            continue
        if seen_zero_row or lead <= last:     # nonzero row below a zero row, or no step to the right
            return False
        last = lead
    return True

print(is_ref([[2, 1, -1, 8], [0, 0.5, 0.5, 1], [0, 0, -1, 1]]))   # True
print(is_ref([[1, 4, 0, 2], [0, 0, 3, 1]]))                       # True
print(is_ref([[1, 2, 3], [0, 0, 0], [0, 0, 4]]))                  # False: zero row not at the bottom
print(is_ref([[1, 2, 3], [0, 5, 1], [0, 3, 2]]))                  # False: two leading entries in column 2

Key takeaways

• REF: zero rows at the bottom, and each leading entry strictly to the right of the one above.
• Every entry below a leading entry is zero, which makes back substitution possible.
• An augmented matrix in REF is inconsistent exactly when a leading entry lies in the right-hand-side column.
• Echelon forms are not unique, but the positions of their leading entries are.

🧪 Practice

1. Which are in REF? (a) $\begin{bmatrix} 0 & 1 & 2 \ 0 & 0 & 3 \end{bmatrix}$ (b) $\begin{bmatrix} 1 & 0 \ 0 & 0 \ 0 & 1 \end{bmatrix}$ (c) $\begin{bmatrix} 3 & 1 & 4 & 1 \ 0 & 0 & 0 & 5 \ 0 & 0 & 0 & 0 \end{bmatrix}$.
2. Solve by back substitution the system whose augmented matrix is $\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \ 0 & 2 & -8 & 8 \ 0 & 0 & 2 & -2 \end{array}\right]$.
3. Give two different echelon forms of $\begin{bmatrix} 2 & 4 \ 1 & 3 \end{bmatrix}$ and check that their leading entries are in the same positions.
4. Interview: Write a function that, given an augmented matrix already in REF, reports "none", "unique", or "infinite" without solving the system. Hint: you only need the column index of each row's leading entry.

Reduced Row Echelon Form

Why it matters. An REF still needs back substitution. The reduced row echelon form shows the solution directly. Unlike the REF, it is unique, so it serves as a fingerprint of a matrix. Two matrices are row equivalent exactly when their reduced forms are equal.

Definition. A matrix is in reduced row echelon form (RREF) if it is in REF and, in addition,

3. every leading entry equals 1, and
4. every leading 1 is the only nonzero entry in its column.

   [ 1  0  *  0  * | * ]       leading entries are 1
   [ 0  1  *  0  * | * ]       zeros above AND below each leading 1
   [ 0  0  0  1  * | * ]       * = any value (only in non-pivot columns
   [ 0  0  0  0  0 | 0 ]           and the right-hand side)

Reading the solution. Each nonzero row of an RREF involves exactly one leading variable, which can be read off directly:

$$\left[\begin{array}{ccc|c} 1 & 0 & 0 & 2 \ 0 & 1 & 0 & 3 \ 0 & 0 & 1 & -1 \end{array}\right] \Rightarrow (x, y, z) = (2, 3, -1); \qquad \left[\begin{array}{ccc|c} 1 & 0 & 2 & 3 \ 0 & 1 & -1 & 4 \end{array}\right] \Rightarrow \begin{cases} x = 3 - 2z \ y = 4 + z \ z \text{ free} \end{cases}$$

Uniqueness theorem. Every matrix is row equivalent to exactly one matrix in reduced row echelon form, written $\operatorname{rref}(A)$.

Idea of the proof. Row operations preserve the solution set of $A\mathbf{x} = \mathbf{0}$. That solution set records every linear relation among the columns of $A$, such as "column 3 equals twice column 1 minus column 2". The RREF is completely determined by these relations. A column is a pivot column exactly when it is not a combination of earlier columns, and the entries of a non-pivot column are the weights that express it in terms of earlier pivot columns. Two RREFs with the same relations are therefore identical.

REF and RREF compared.

	REF	RREF
Unique for a given matrix	No	Yes
Leading entries	Any nonzero value	Exactly 1
Entries above leading entries	Arbitrary	Zero
Solving requires	Back substitution	Reading off
Cost to reach (square, $n \times n$)	$\approx \tfrac23 n^3$ operations	$\approx n^3$ operations
Typical use	Numerical solving	Hand work, theory, describing solution sets

import sympy as sp

M = sp.Matrix([[2, 1, -1, 8], [-3, -1, 2, -11], [-2, 1, 2, -3]])
R, pivots = M.rref()                  # exact arithmetic over the rationals
print(R)                              # Matrix([[1, 0, 0, 2], [0, 1, 0, 3], [0, 0, 1, -1]])
print(pivots)                         # (0, 1, 2)

# Different-looking matrices with the same RREF are row equivalent
P = sp.Matrix([[1, 2, 3], [2, 4, 7]])
Q = sp.Matrix([[3, 6, 10], [1, 2, 2]])
print(P.rref()[0])                    # Matrix([[1, 2, 0], [0, 0, 1]])
print(P.rref()[0] == Q.rref()[0])     # True

Key takeaways

• RREF adds two conditions to REF: every leading entry is 1, and it is the only nonzero entry in its column.
• In RREF, each nonzero row solves for one leading variable in terms of the others.
• The RREF of a matrix is unique, so two matrices are row equivalent exactly when they have the same RREF.
• RREF costs more to compute than REF; it is best for hand work, theory, and describing solution sets.

🧪 Practice

1. Reduce $\begin{bmatrix} 1 & 2 & 1 \ 2 & 5 & 3 \end{bmatrix}$ to RREF.
2. Is $\begin{bmatrix} 1 & 3 & 0 & -1 \ 0 & 0 & 1 & 2 \ 0 & 0 & 0 & 0 \end{bmatrix}$ in RREF? Interpreting it as an augmented matrix, write the general solution.
3. Decide whether $\begin{bmatrix} 1 & 1 \ 2 & 3 \end{bmatrix}$ and $\begin{bmatrix} 0 & 1 \ 5 & 0 \end{bmatrix}$ are row equivalent by computing both RREFs.
4. Interview: How many distinct $2 \times 2$ matrices are in reduced row echelon form, if you count each family with free parameters as one pattern? List the patterns. Hint: classify by the number of pivots (0, 1, or 2) and by where they can sit.

Gaussian Elimination Algorithm

Why it matters. Gaussian elimination is a finite, mechanical procedure that decides consistency and finds every solution of every linear system. With refinements, it is what numerical libraries such as LAPACK run when you call a linear solver.

The forward phase (to REF). Work on the augmented matrix, column by column, from left to right:

1. Start with the current row $r$ at the top and the current column $c$ at the far left.
2. In column $c$, look for a nonzero entry in row $r$ or below. If there is none, move to the next column and repeat this step.
3. Interchange rows to bring that entry into row $r$. This entry is the pivot.
4. Use replacements $R_i \leftarrow R_i - \frac{a_{ic}}{a_{rc}}R_r$ to create zeros below the pivot. The ratio $\frac{a_{ic}}{a_{rc}}$ is the multiplier.
5. Move to row $r + 1$ and column $c + 1$, and repeat until you run out of rows or columns.

The backward phase. If a row $[,0 \cdots 0 \mid c,]$ with $c \ne 0$ appears, stop: the system is inconsistent. Otherwise, solve by back substitution, or continue to RREF (Gauss-Jordan, next topic).

Example. The full forward phase on the running example:

 [  2   1  -1 |   8 ]        [ 2   1    -1   |  8 ]        [ 2   1    -1   |  8 ]
 [ -3  -1   2 | -11 ]  --->  [ 0   1/2  1/2  |  1 ]  --->  [ 0   1/2  1/2  |  1 ]
 [ -2   1   2 |  -3 ]        [ 0   2    1    |  5 ]        [ 0   0    -1   |  1 ]

   pivot 2 in column 1:        pivot 1/2 in column 2:
   R2 <- R2 + (3/2) R1         R3 <- R3 - 4 R2
   R3 <- R3 + R1

Back substitution gives $z = -1$, $y = 3$, $x = 2$.

Partial pivoting. In exact arithmetic, any nonzero pivot works. In floating point, a tiny pivot produces huge multipliers that swamp the other entries in rounding error. Partial pivoting chooses, in step 2, the candidate with the largest absolute value. Consider

$$\left[\begin{array}{cc|c} 10^{-20} & 1 & 1 \ 1 & 1 & 2 \end{array}\right], \qquad \text{exact solution} \approx (1, 1).$$

Without a row swap, the multiplier is $10^{20}$. Row 2 becomes $[,0 \ \ 1 - 10^{20} \mid 2 - 10^{20},]$, which rounds to $[,0 \ \ {-10^{20}} \mid -10^{20},]$. That gives $y = 1$ and then $x = (1 - 1)/10^{-20} = 0$, which is completely wrong. Swapping the rows first gives the correct answer. Chapter 12 studies stability in detail.

Cost. For an $n \times n$ system, the forward phase takes about $\tfrac23 n^3$ floating-point operations and back substitution about $n^2$:

$n$	Approximate operations	Time at $10^9$ operations per second
10	$7 \times 10^2$	under a microsecond
1,000	$7 \times 10^8$	under a second
100,000	$7 \times 10^{14}$	about 8 days

The cubic growth is why large sparse systems call for iterative methods (Section 12.2).

import numpy as np

def gaussian_elimination(A, b, pivoting=True):
    """Solve a square system Ax = b by forward elimination and back substitution."""
    M = np.column_stack([A, b]).astype(float)       # work on a copy of [A | b]
    n = len(b)
    for k in range(n):                              # k = current pivot row and column
        if pivoting:
            p = k + np.argmax(np.abs(M[k:, k]))     # largest candidate in column k, at or below row k
            M[[k, p]] = M[[p, k]]                   # interchange rows k and p
        if M[k, k] == 0:
            raise ValueError("zero pivot: no unique solution")
        for i in range(k + 1, n):
            m = M[i, k] / M[k, k]                   # multiplier
            M[i, k:] -= m * M[k, k:]                # replacement R_i <- R_i - m R_k
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):                  # back substitution
        x[i] = (M[i, n] - M[i, i+1:n] @ x[i+1:]) / M[i, i]
    return x

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
b = np.array([8., -11, -3])
print(gaussian_elimination(A, b))                   # approximately [ 2.  3. -1.]

# A tiny pivot destroys accuracy unless the rows are swapped first
A2 = np.array([[1e-20, 1.], [1., 1.]])
b2 = np.array([1., 2.])
print(gaussian_elimination(A2, b2, pivoting=False)) # [0. 1.]: x is completely wrong
print(gaussian_elimination(A2, b2, pivoting=True))  # [1. 1.]

Key takeaways

• The forward phase moves left to right, choosing a pivot in each column and creating zeros below it, until the matrix is in REF.
• A row $[,0 \cdots 0 \mid c,]$ with $c \ne 0$ means inconsistency; otherwise back substitution finishes the job.
• In floating point, partial pivoting (choosing the largest available pivot) is essential for accuracy.
• Elimination costs about $\tfrac23 n^3$ operations for an $n \times n$ system, which limits dense direct methods to moderate $n$.

🧪 Practice

1. Solve by Gaussian elimination: $x + y + z = 2$, $\ 2x + 3y + z = 3$, $\ x - y - 2z = -6$.
2. Apply Gaussian elimination to $x + 2y - z = 1$, $\ 2x + 4y + z = 5$, $\ 3x + 6y = 6$, and describe what happens at the second column.
3. Count the multiplications and divisions performed in the forward phase for a general $3 \times 3$ system (augmented $3 \times 4$). Compare with $\tfrac23 n^3$ for $n = 3$.
4. Interview: Why does partial pivoting choose the largest pivot rather than any nonzero one? Explain what goes wrong numerically when a pivot is tiny. Hint: look at the size of the multipliers and at what happens when a huge number is added to a moderate one in floating point.

Gauss-Jordan Elimination

Why it matters. Gauss-Jordan elimination continues past REF all the way to RREF. It needs no separate back substitution, handles every shape of system including those with free variables, and is the standard method for computing inverses by hand ($[,A \mid I,] \to [,I \mid A^{-1},]$, Section 4.3).

Algorithm. First run the forward phase of Gaussian elimination. Then, in a backward phase, start with the rightmost pivot and move up and to the left:

1. Scale the pivot's row so that the pivot becomes 1.
2. Use replacements to create zeros above the pivot.

A common variant scales each pivot to 1 and clears its column both above and below in a single forward pass. The result is the same, because the RREF is unique.

Example. Continuing the running example from its REF:

 [ 2   1    -1   |  8 ]   R3 <- -R3         [ 2  1  -1 |  8 ]   R2 <- R2 - R3   [ 2  1  0 |  7 ]
 [ 0   1/2  1/2  |  1 ]   R2 <- 2 R2        [ 0  1   1 |  2 ]   R1 <- R1 + R3   [ 0  1  0 |  3 ]
 [ 0   0    -1   |  1 ]  ------------->     [ 0  0   1 | -1 ]  ------------->   [ 0  0  1 | -1 ]

                          R1 <- R1 - R2     [ 2  0  0 |  4 ]   R1 <- R1 / 2    [ 1  0  0 |  2 ]
                         ------------->     [ 0  1  0 |  3 ]  ------------->   [ 0  1  0 |  3 ]
                                            [ 0  0  1 | -1 ]                   [ 0  0  1 | -1 ]

The solution $(2, 3, -1)$ appears in the last column.

Gaussian versus Gauss-Jordan.

	Gaussian elimination + back substitution	Gauss-Jordan elimination
End form	REF	RREF
Operations ($n \times n$)	$\approx \tfrac23 n^3$	$\approx n^3$
Solution	Computed by back substitution	Read off the last column
Free variables	Handled during back substitution	Visible directly
Typical use	Numerical software (via LU, Section 4.4)	Hand computation, inverses, theory

from fractions import Fraction as F

def rref(M):
    """Gauss-Jordan elimination with exact fractions. Returns (RREF, list of pivot columns)."""
    A = [[F(v) for v in row] for row in M]
    rows, cols = len(A), len(A[0])
    pivots, r = [], 0
    for c in range(cols):
        p = next((i for i in range(r, rows) if A[i][c] != 0), None)
        if p is None:
            continue                                   # no pivot in this column: move right
        A[r], A[p] = A[p], A[r]                        # interchange
        A[r] = [v / A[r][c] for v in A[r]]             # scale the pivot to 1
        for i in range(rows):                          # clear the column above and below the pivot
            factor = A[i][c]
            if i != r and factor != 0:
                A[i] = [a - factor * ar for a, ar in zip(A[i], A[r])]
        pivots.append(c)
        r += 1
        if r == rows:
            break
    return A, pivots

def show(A):
    for row in A:
        print("  ".join(f"{str(v):>4}" for v in row))

R, piv = rref([[2, 1, -1, 8], [-3, -1, 2, -11], [-2, 1, 2, -3]])
show(R)            # last column: 2, 3, -1
print(piv)         # [0, 1, 2]

# A system with more unknowns than pivots: 3 equations, 4 unknowns (augmented 3 x 5)
R, piv = rref([[1, 2, -1, 3, 4], [2, 4, 0, 2, 6], [1, 2, 1, -1, 2]])
show(R)            # rows: [1 2 0 1 | 3], [0 0 1 -2 | -1], [0 0 0 0 | 0]
print(piv)         # [0, 2]

Key takeaways

• Gauss-Jordan elimination reduces all the way to RREF: each pivot is scaled to 1, and its column is cleared above and below.
• The solution, or the general solution with free variables, is read directly from the RREF.
• It costs about $n^3$ operations versus $\tfrac23 n^3$ for Gaussian elimination with back substitution, so numerical code prefers the latter.
• Applied to $[,A \mid I,]$, Gauss-Jordan computes the inverse of $A$ (Section 4.3).

🧪 Practice

1. Use Gauss-Jordan elimination to solve $x + 2y = 4$, $\ 3x - y = 5$.
2. Reduce $\left[\begin{array}{ccc|c} 1 & 1 & 1 & 3 \ 1 & 2 & 3 & 6 \ 2 & 3 & 4 & 9 \end{array}\right]$ to RREF and write the general solution.
3. Apply Gauss-Jordan to $[,A \mid I,]$ for $A = \begin{bmatrix} 2 & 1 \ 5 & 3 \end{bmatrix}$, and verify that the right half is a matrix $B$ with $AB = I$.
4. Interview: Library routines that "solve $A\mathbf{x} = \mathbf{b}$" almost never compute the RREF or the inverse of $A$. Why not? Hint: compare operation counts, and think about how rounding errors accumulate in the extra work.

Pivot Positions and Pivot Columns

Why it matters. Pivot positions hold the structural information about a matrix. Their number is the rank. Their positions decide which variables are basic and which are free, whether a system is consistent, whether solutions are unique, and which columns form a basis of the column space (Section 6.6).

Definitions.

• A pivot position of $A$ is a location $(i, j)$ that holds a leading 1 in $\operatorname{rref}(A)$.
• A pivot column is a column of $A$ that contains a pivot position.
• A pivot is the nonzero number placed in a pivot position during elimination. Its value depends on the row operations used, but its position does not.

Pivot positions are well defined because the RREF is unique. Any REF has its leading entries in the same positions, since the backward phase from REF to RREF never moves a leading entry.

Example. For the coefficient matrix from the Gauss-Jordan example,

       column:  1    2    3    4
  A        = [  1    2   -1    3 ]          rref(A) = [ 1    2    0    1 ]   <- pivot in column 1
             [  2    4    0    2 ]                    [ 0    0    1   -2 ]   <- pivot in column 3
             [  1    2    1   -1 ]                    [ 0    0    0    0 ]   <- no pivot in row 3
                                                        ^         ^
                                       pivot columns of A: 1 and 3;   non-pivot columns: 2 and 4

The pivot columns are columns of the original matrix $A$. Columns of $\operatorname{rref}(A)$ have different entries. Only the column positions carry over.

What pivots tell you. For a system $A\mathbf{x} = \mathbf{b}$ with $A$ of size $m \times n$:

Question	Answer in terms of pivots
Is $A\mathbf{x} = \mathbf{b}$ consistent?	Yes, if the last column of $[,A \mid \mathbf{b},]$ is not a pivot column
If consistent, is the solution unique?	Yes, if every column of $A$ is a pivot column
Is $A\mathbf{x} = \mathbf{b}$ consistent for every $\mathbf{b}$?	Yes, if every row of $A$ has a pivot
What is the rank of $A$?	The number of pivot positions

Counting consequences. Each row and each column holds at most one pivot, so the number of pivots is at most $\min(m, n)$. It follows that:

• if $m < n$, some column of $A$ has no pivot, so a consistent system has infinitely many solutions;
• if $m > n$, some row of $A$ has no pivot, so some right-hand side $\mathbf{b}$ gives an inconsistent system.

import sympy as sp

A = sp.Matrix([[1, 2, -1, 3], [2, 4, 0, 2], [1, 2, 1, -1]])
b = sp.Matrix([4, 6, 2])
m, n = A.shape

_, piv_A = A.rref()
_, piv_Ab = A.row_join(b).rref()                    # pivots of the augmented matrix [A | b]
print(piv_A)                                        # (0, 2): pivot columns 1 and 3
print("consistent:", n not in piv_Ab)               # True: no pivot in the augmented column
print("unique:", len(piv_A) == n)                   # False: columns 2 and 4 have no pivot
print("solvable for every b:", len(piv_A) == m)     # False: row 3 has no pivot
print("pivot columns of A:", [list(A.col(j)) for j in piv_A])   # taken from the original A

Key takeaways

• Pivot positions are the locations of the leading 1s in $\operatorname{rref}(A)$; they are the same in every echelon form.
• Pivot columns are columns of the original $A$ that contain a pivot position.
• Consistency, uniqueness, solvability for every $\mathbf{b}$, and rank can all be read from the pivot positions.
• There are at most $\min(m, n)$ pivots, so wide systems are never uniquely solvable and tall systems are not solvable for every $\mathbf{b}$.

🧪 Practice

1. Find the pivot columns of $\begin{bmatrix} 0 & 3 & -6 & 6 \ 3 & -7 & 8 & -5 \ 3 & -9 & 12 & -9 \end{bmatrix}$.
2. A $5 \times 7$ coefficient matrix has 5 pivots. Is $A\mathbf{x} = \mathbf{b}$ consistent for every $\mathbf{b}$? Is the solution unique? Explain using rows and columns.
3. A $6 \times 4$ coefficient matrix has 4 pivots. What can you say about the number of solutions of $A\mathbf{x} = \mathbf{b}$ for different $\mathbf{b}$?
4. Interview: Two engineers reduce the same matrix with different row operations and get different echelon forms. One claims the pivot columns therefore depend on the method. Refute the claim. Hint: relate every echelon form to the unique RREF.

<a id="34-structure-of-solution-sets"></a>

3.4 Structure of Solution Sets

Once a system is in reduced echelon form, its solution set can be written down explicitly; this subchapter shows how, and proves that every solution set is one particular solution plus all solutions of the associated homogeneous system.

Throughout this subchapter, the running example is the $3 \times 4$ system reduced in Section 3.3:

$$ \left[\begin{array}{cccc|c} 1 & 2 & -1 & 3 & 4 \ 2 & 4 & 0 & 2 & 6 \ 1 & 2 & 1 & -1 & 2 \end{array}\right] \ \xrightarrow{\ \text{rref}\ }
\left[\begin{array}{cccc|c} 1 & 2 & 0 & 1 & 3 \ 0 & 0 & 1 & -2 & -1 \ 0 & 0 & 0 & 0 & 0 \end{array}\right]. $$

Free Variables and Basic Variables

Why it matters. When a solution is not unique, you need a way to describe all solutions. Splitting the variables into those you may choose (free) and those that are then forced (basic) gives that description. The number of free variables is the number of degrees of freedom in the solution set.

Definitions. For a consistent system with coefficient matrix $A$:

• a basic variable (or leading variable) corresponds to a pivot column of $A$;
• a free variable corresponds to a non-pivot column of $A$.

In the RREF, each nonzero row contains exactly one basic variable. Solve that row for its basic variable in terms of the free ones.

Example. The RREF above reads

$$ \begin{aligned} x_1 + 2x_2 \phantom{{}+x_3} + x_4 &= 3 \ x_3 - 2x_4 &= -1 \end{aligned} \quad\Longrightarrow\quad \begin{cases} x_1 = 3 - 2x_2 - x_4 \ x_2 \text{ free} \ x_3 = -1 + 2x_4 \ x_4 \text{ free} \end{cases} $$

The pivot columns are 1 and 3, so $x_1$ and $x_3$ are basic and $x_2$ and $x_4$ are free. Each choice of $(x_2, x_4)$ produces exactly one solution. Think of a mechanism with two control knobs: once you set the knobs $x_2$ and $x_4$, every other part's position is determined.

Counting.

Quantity	Equals
Number of basic variables	Number of pivots, $r$
Number of free variables	$n - r$
Dimension of the solution set (if consistent)	$n - r$

Two cautions.

• The labels depend on column order; the count does not. For $x + y = 1$ with order $(x, y)$, $y$ is free. With order $(y, x)$, $x$ is free. Either way there is one free variable.
• Free variables only matter for consistent systems. If $[,A \mid \mathbf{b},]$ has a pivot in its last column, there are no solutions, whatever $A$'s free columns suggest.

import sympy as sp

aug = sp.Matrix([[1, 2, -1, 3, 4], [2, 4, 0, 2, 6], [1, 2, 1, -1, 2]])
R, pivots = aug.rref()
n = aug.shape[1] - 1                                  # number of unknowns
xs = sp.symbols(f"x1:{n + 1}")                        # x1, x2, x3, x4

if n in pivots:
    print("inconsistent")
else:
    free = [xs[j] for j in range(n) if j not in pivots]
    print("basic:", [xs[p] for p in pivots], "free:", free)   # basic: [x1, x3] free: [x2, x4]
    for i, p in enumerate(pivots):                    # row i of the RREF determines basic variable x_p
        rhs = R[i, n] - sum(R[i, j] * xs[j] for j in range(n) if j != p)
        print(f"{xs[p]} = {rhs}")
# x1 = -2*x2 - x4 + 3
# x3 = 2*x4 - 1

Key takeaways

• Basic variables belong to pivot columns; free variables belong to non-pivot columns.
• Each nonzero row of the RREF expresses one basic variable in terms of the free variables.
• A consistent system with $n$ unknowns and $r$ pivots has $n - r$ free variables, and that is the dimension of its solution set.
• Which variables are free depends on the column order, but how many are free does not.

🧪 Practice

1. Identify the basic and free variables and write the general solution for $\left[\begin{array}{cccc|c} 1 & 0 & -3 & 0 & 2 \ 0 & 1 & 4 & 0 & -1 \ 0 & 0 & 0 & 1 & 5 \end{array}\right]$.
2. Solve $x_1 - x_2 + 2x_3 = 1$, $\ 2x_1 - 2x_2 + 5x_3 = 4$ and identify the free variables.
3. Reorder the unknowns of the running example as $x_4, x_3, x_2, x_1$ and reduce again. Which variables are free now?
4. Interview: A system has 8 unknowns and its coefficient matrix has 5 pivots. A teammate says "the system has 3 degrees of freedom." Is that always true? Hint: the claim assumes something about the augmented column.

Parametric Vector Form of Solutions

Why it matters. Listing the basic variables in terms of the free ones describes the solutions, but it hides their geometry. Parametric vector form rewrites the same description as a point plus a span of direction vectors. That is exactly the vector form of lines and planes from Section 2.5, now in any dimension. It also makes generating, checking, and reasoning about solutions straightforward.

Procedure.

1. Write the solution vector $\mathbf{x}$, replacing each basic variable by its expression in the free variables.
2. Rename the free variables as parameters, for example $x_2 = s$ and $x_4 = t$.
3. Split the vector into a constant part plus one vector per parameter.

Example. From $x_1 = 3 - 2s - t$, $x_2 = s$, $x_3 = -1 + 2t$, $x_4 = t$:

        [ x1 ]   [ 3 - 2s - t ]   [  3 ]     [ -2 ]     [ -1 ]
    x = [ x2 ] = [     s      ] = [  0 ] + s [  1 ] + t [  0 ]
        [ x3 ]   [  -1 + 2t   ]   [ -1 ]     [  0 ]     [  2 ]
        [ x4 ]   [     t      ]   [  0 ]     [  0 ]     [  1 ]
                                     p          v1         v2

$$\mathbf{x} = \mathbf{p} + s\mathbf{v}_1 + t\mathbf{v}_2, \qquad s, t \in \mathbb{R}.$$

The solution set is a 2-dimensional plane in $\mathbb{R}^4$ through $\mathbf{p} = (3, 0, -1, 0)$, spanned by $\mathbf{v}_1$ and $\mathbf{v}_2$.

Structure of the direction vectors. Each $\mathbf{v}_j$ has a 1 in the position of its own free variable and 0 in the positions of the other free variables. For example, $\mathbf{v}_1$ has entries $(x_2, x_4) = (1, 0)$ and $\mathbf{v}_2$ has $(0, 1)$. So no direction vector is a combination of the others, and different parameter values always give different solutions.

Checking the result. Two checks catch nearly all mistakes: $A\mathbf{p} = \mathbf{b}$, and $A\mathbf{v}_j = \mathbf{0}$ for each $j$. The second check anticipates the next topic: the direction vectors solve the homogeneous system.

Representations are not unique. Another particular point, such as $\mathbf{p} + \mathbf{v}_1$, or other spanning directions, such as $\mathbf{v}_1 + \mathbf{v}_2$ and $\mathbf{v}_2$, describe the same set. This is the same freedom you have when writing the equation of a line.

import numpy as np

A = np.array([[1., 2, -1, 3], [2, 4, 0, 2], [1, 2, 1, -1]])
b = np.array([4., 6, 2])
p  = np.array([3., 0, -1, 0])
v1 = np.array([-2., 1, 0, 0])
v2 = np.array([-1., 0, 2, 1])

print(A @ p, A @ v1, A @ v2)          # [4. 6. 2.] [0. 0. 0.] [0. 0. 0.]

rng = np.random.default_rng(0)
for s, t in rng.normal(size=(3, 2)):  # any parameter values give a solution
    x = p + s * v1 + t * v2
    print(np.allclose(A @ x, b))      # True, True, True

Key takeaways

• Parametric vector form writes the general solution as $\mathbf{x} = \mathbf{p} + t_1\mathbf{v}_1 + \cdots + t_k\mathbf{v}_k$, one parameter per free variable.
• The solution set is a $k$-dimensional flat: a point, line, plane, or higher analog, as in Section 2.5.
• Each direction vector has a 1 in its own free-variable slot and 0 in the other free-variable slots.
• Verify with $A\mathbf{p} = \mathbf{b}$ and $A\mathbf{v}_j = \mathbf{0}$; other valid parametrizations of the same set exist.

🧪 Practice

1. Write the solution of $x_1 + 3x_2 - x_3 = 2$ (one equation, three unknowns) in parametric vector form. What shape is it?
2. Write the general solution of $\left[\begin{array}{cccc|c} 1 & -2 & 0 & 3 & 1 \ 0 & 0 & 1 & -1 & 4 \end{array}\right]$ in parametric vector form and verify both checks.
3. Show that $\mathbf{x} = (1, 1, -1, 0) + s(-2, 1, 0, 0) + t(-3, 1, 2, 1)$ describes the same solution set as the running example.
4. Interview: A recommendation system needs any vector $\mathbf{x}$ with $A\mathbf{x} = \mathbf{b}$, and among all such vectors, the one with the smallest norm. Given the parametric vector form, how would you find it? Hint: the squared norm $\lVert \mathbf{p} + \sum t_j\mathbf{v}_j \rVert^2$ is a quadratic function of the parameters; alternatively, look for the solution orthogonal to every $\mathbf{v}_j$.

Homogeneous Systems

Why it matters. The system $A\mathbf{x} = \mathbf{0}$ is the version of $A\mathbf{x} = \mathbf{b}$ with every constant set to zero. It always has a solution, and its solution set measures how far solutions of every system with matrix $A$ are from unique. Homogeneous systems also appear directly: in balancing chemical equations (Section 3.5), in testing linear independence (Section 6.3), and in finding eigenvectors from $(A - \lambda I)\mathbf{x} = \mathbf{0}$ (Chapter 8).

Definition. A system is homogeneous if it can be written $A\mathbf{x} = \mathbf{0}$. The zero vector is always a solution, called the trivial solution. A nontrivial solution is any solution $\mathbf{x} \ne \mathbf{0}$.

When nontrivial solutions exist. Row operations never change a zero column, so the augmented column stays zero and the system is always consistent. Then:

$$A\mathbf{x} = \mathbf{0} \text{ has a nontrivial solution} \iff \text{the system has at least one free variable}.$$

Theorem. A homogeneous system with more unknowns than equations ($n > m$) has a nontrivial solution.

Proof. There are at most $m$ pivots, and $m < n$, so some column has no pivot. That column gives a free variable. Setting it to 1 produces a nonzero solution. $\blacksquare$

Closure. If $\mathbf{u}$ and $\mathbf{v}$ solve $A\mathbf{x} = \mathbf{0}$, then so do $\mathbf{u} + \mathbf{v}$ and $c\mathbf{u}$:

$$A(\mathbf{u} + \mathbf{v}) = A\mathbf{u} + A\mathbf{v} = \mathbf{0} + \mathbf{0} = \mathbf{0}, \qquad A(c\mathbf{u}) = c,A\mathbf{u} = \mathbf{0}.$$

So the solution set is closed under addition and scalar multiplication and contains $\mathbf{0}$. It is a subspace (Section 6.2), called the null space of $A$, written $\operatorname{Nul} A$. In parametric vector form it has no constant term:

$$\operatorname{Nul} A = {, t_1\mathbf{v}_1 + \cdots + t_k\mathbf{v}_k ,} = \operatorname{Span}{\mathbf{v}_1, \dots, \mathbf{v}_k}.$$

Example. For the running example's matrix, $\operatorname{Nul} A = \operatorname{Span}{(-2, 1, 0, 0),\ (-1, 0, 2, 1)}$, a plane through the origin in $\mathbb{R}^4$. For the smaller system $x + y + z = 0$, $x - y = 0$, the solutions satisfy $x = y$ and $z = -2x$, so $\operatorname{Nul} = \operatorname{Span}{(1, 1, -2)}$. Two planes through the origin meet in a line through the origin.

    homogeneous solution sets always pass through the origin:

        z
        |     .  line Span{(1, 1, -2)}
        |   .
        | .
        O-------------- y          (planes through O meet in a line through O)
       / .
      /    .
     x       .

Property	$A\mathbf{x} = \mathbf{0}$	$A\mathbf{x} = \mathbf{b}$, $\ \mathbf{b} \ne \mathbf{0}$
Always consistent	Yes ($\mathbf{x} = \mathbf{0}$)	No
Contains the origin	Yes	No
Closed under $+$ and scaling	Yes: a subspace	No
Parametric form	$t_1\mathbf{v}_1 + \cdots + t_k\mathbf{v}_k$	$\mathbf{p} + t_1\mathbf{v}_1 + \cdots + t_k\mathbf{v}_k$

import sympy as sp

A = sp.Matrix([[1, 2, -1, 3], [2, 4, 0, 2], [1, 2, 1, -1]])
basis = A.nullspace()                          # exact basis of the solutions of Ax = 0
for v in basis:
    print(v.T, (A * v).T)                      # each basis vector is sent to the zero vector
# Matrix([[-2, 1, 0, 0]]) Matrix([[0, 0, 0]])
# Matrix([[-1, 0, 2, 1]]) Matrix([[0, 0, 0]])

u = 3 * basis[0] - 5 * basis[1]                # closure: combinations of solutions are solutions
print((A * u).T)                               # Matrix([[0, 0, 0]])

B = sp.Matrix([[1, 1, 1], [1, -1, 0]])         # 2 equations, 3 unknowns: a nontrivial solution must exist
print(B.nullspace())                           # [Matrix([[-1/2], [-1/2], [1]])], a multiple of (1, 1, -2)

Key takeaways

• $A\mathbf{x} = \mathbf{0}$ is always consistent; $\mathbf{x} = \mathbf{0}$ is the trivial solution.
• Nontrivial solutions exist exactly when there is a free variable, which is guaranteed whenever $n > m$.
• Sums and scalar multiples of homogeneous solutions are solutions, so the solution set is a subspace: the null space $\operatorname{Nul} A$.
• Since the zero right-hand side never changes, you can row reduce $A$ alone.

🧪 Practice

1. Find a basis for the solution set of $x_1 - 3x_2 + 2x_3 = 0$ (one equation in $\mathbb{R}^3$) and describe it geometrically.
2. Determine whether $\begin{bmatrix} 1 & 2 & 3 \ 4 & 5 & 6 \ 7 & 8 & 9 \end{bmatrix}\mathbf{x} = \mathbf{0}$ has a nontrivial solution. If so, find one.
3. Prove that the solution set of $A\mathbf{x} = \mathbf{b}$ with $\mathbf{b} \ne \mathbf{0}$ is never a subspace.
4. Interview: In a data pipeline, the columns of a feature matrix $X$ are features and the rows are samples. You discover that $X\mathbf{w} = \mathbf{0}$ has a nontrivial solution $\mathbf{w}$. What does that say about the features, and why is it a problem for fitting a linear model? Hint: read $X\mathbf{w} = \mathbf{0}$ as a statement about a combination of columns, and ask whether the fitted weights can still be unique.

Nonhomogeneous Systems and Particular Solutions

Why it matters. The structure theorem below says that every solution of $A\mathbf{x} = \mathbf{b}$ is one particular solution plus a solution of $A\mathbf{x} = \mathbf{0}$. It splits every problem into two independent parts: find any single solution, then describe the ambiguity, which depends only on $A$. The same principle, "general solution = particular + homogeneous", governs linear differential equations and linear recurrences.

Theorem. Suppose $A\mathbf{x} = \mathbf{b}$ is consistent, and let $\mathbf{p}$ be any one solution, called a particular solution. Then the solution set is

$${, \mathbf{x} : A\mathbf{x} = \mathbf{b} ,} = {, \mathbf{p} + \mathbf{v}_h : A\mathbf{v}_h = \mathbf{0} ,} = \mathbf{p} + \operatorname{Nul} A.$$

Proof. If $A\mathbf{v}_h = \mathbf{0}$, then $A(\mathbf{p} + \mathbf{v}_h) = \mathbf{b} + \mathbf{0} = \mathbf{b}$. Conversely, if $A\mathbf{y} = \mathbf{b}$, then $A(\mathbf{y} - \mathbf{p}) = \mathbf{b} - \mathbf{b} = \mathbf{0}$, so $\mathbf{y} = \mathbf{p} + (\mathbf{y} - \mathbf{p})$ has the required form. $\blacksquare$

Geometry. The solution set of $A\mathbf{x} = \mathbf{b}$ is the null space shifted by $\mathbf{p}$. It is parallel to $\operatorname{Nul} A$ but generally does not pass through the origin.

                Nul A                          p + Nul A
                  /                               /
                 /                               /
                /                               /
               /             p                 /
              O  - - - - - - - - - - - - - - -*
             /                                /
            /                                /
    (solutions of Ax = 0)           (solutions of Ax = b)

An analogy: think of $A\mathbf{x}$ as a sensor reading. Directions in $\operatorname{Nul} A$ are invisible to the sensor, because moving along them does not change the reading. All inputs that produce the reading $\mathbf{b}$ are therefore one such input plus any invisible movement.

Consequences.

Fact	Reason
A consistent system $A\mathbf{x} = \mathbf{b}$ has a unique solution exactly when $A\mathbf{x} = \mathbf{0}$ has only the trivial solution	The solution set is $\mathbf{p} + \operatorname{Nul} A$
Uniqueness does not depend on $\mathbf{b}$	$\operatorname{Nul} A$ depends only on $A$
Any particular solution works	Two particular solutions differ by an element of $\operatorname{Nul} A$
Superposition: if $A\mathbf{p}_1 = \mathbf{b}_1$ and $A\mathbf{p}_2 = \mathbf{b}_2$, then $A(\mathbf{p}_1 + \mathbf{p}_2) = \mathbf{b}_1 + \mathbf{b}_2$	Linearity of $A$

Example. In the running example, $\mathbf{p} = (3, 0, -1, 0)$ came from setting the free variables to zero. Setting $s = 1$, $t = 0$ instead gives $\mathbf{p}' = (1, 1, -1, 0)$. Their difference $\mathbf{p}' - \mathbf{p} = (-2, 1, 0, 0) = \mathbf{v}_1$ lies in $\operatorname{Nul} A$, as the theorem predicts.

import numpy as np
import scipy.linalg as la

A = np.array([[1., 2, -1, 3], [2, 4, 0, 2], [1, 2, 1, -1]])
b = np.array([4., 6, 2])

p, *_ = np.linalg.lstsq(A, b, rcond=None)    # one particular solution (here, the one of smallest norm)
N = la.null_space(A)                         # columns: an orthonormal basis of Nul A, shape (4, 2)
print(np.allclose(A @ p, b), N.shape)        # True (4, 2)

c = np.array([2.5, -1.0])                    # any coefficients give another solution p + N c
print(np.allclose(A @ (p + N @ c), b))       # True

p_hand = np.array([3., 0, -1, 0])            # the particular solution found by hand
print(np.allclose(A @ (p_hand - p), 0))      # True: two particular solutions differ by a null-space vector

Key takeaways

• If $\mathbf{p}$ is one solution of $A\mathbf{x} = \mathbf{b}$, then every solution is $\mathbf{p} + \mathbf{v}_h$ with $\mathbf{v}_h \in \operatorname{Nul} A$.
• Geometrically, the solution set is the null space translated by $\mathbf{p}$; it contains the origin only when $\mathbf{b} = \mathbf{0}$.
• A consistent system has a unique solution exactly when $\operatorname{Nul} A = {\mathbf{0}}$, regardless of $\mathbf{b}$.
• Any particular solution can be used, and solutions for different right-hand sides add (superposition).

🧪 Practice

1. Given that $(1, 2, 0)$ solves $A\mathbf{x} = \mathbf{b}$ and $\operatorname{Nul} A = \operatorname{Span}{(1, -1, 1)}$, describe all solutions and give two more.
2. Find a particular solution of $x + y + z = 3$, $\ x - y = 1$, find $\operatorname{Nul}$ of the coefficient matrix, and combine them into the general solution.
3. If $A\mathbf{p}_1 = (1, 0)$ and $A\mathbf{p}_2 = (0, 1)$, find a solution of $A\mathbf{x} = (3, -2)$.
4. Interview: A linear system models a mechanical structure, and its matrix has a nontrivial null space. What does a nonzero vector in that null space represent physically, and why do engineers worry about it? Hint: think of a displacement that produces no change in the forces: a motion the structure cannot resist.

Rank and Number of Solutions

Why it matters. One number, the rank, measures how many genuinely independent equations a system contains. Comparing $\operatorname{rank} A$, $\operatorname{rank},[,A \mid \mathbf{b},]$, and $n$ classifies every linear system without solving it.

Definition. The rank of a matrix $A$ is the number of pivot positions in $A$, which equals the number of nonzero rows in any echelon form of $A$. It is well defined because pivot positions do not depend on the row operations used. Two interpretations, developed in Chapter 6:

• the number of independent equations, after redundant ones are discarded;
• the dimension of the span of the columns of $A$, the column space.

From the previous topics, the number of free variables is $n - \operatorname{rank} A$, which is the dimension of $\operatorname{Nul} A$. The identity

$$\operatorname{rank} A + \dim \operatorname{Nul} A = n$$

is the rank-nullity theorem (Sections 6.6 and 7.3).

Rouché-Capelli theorem. For an $m \times n$ system, let $r = \operatorname{rank} A$ and $r' = \operatorname{rank},[,A \mid \mathbf{b},]$. Adding one column can add at most one pivot, so $r' = r$ or $r' = r + 1$.

Condition	Outcome
$r' = r + 1$	No solution: the augmented column is a pivot column
$r' = r = n$	Exactly one solution
$r' = r < n$	Infinitely many solutions, with $n - r$ free parameters

              compute  r = rank A   and   r' = rank [A | b]
                                  |
                       r' > r ?  -+-- yes -->  no solution
                                  |
                                  no
                                  |
                       r = n ?   -+-- yes -->  exactly one solution
                                  |
                                  no
                                  |
                                  +--------->  infinitely many: n - r free parameters

Rank by shape. Because $r \le \min(m, n)$:

Shape	What rank allows	Consequence
Wide, $m < n$	$r \le m < n$	Never a unique solution; $A\mathbf{x} = \mathbf{0}$ has nontrivial solutions
Square, $m = n$	$r = n$ is possible	If $r = n$: exactly one solution for every $\mathbf{b}$ (invertible, Section 4.3)
Tall, $m > n$	$r \le n < m$	Some $\mathbf{b}$ give no solution; consistent systems are unique when $r = n$
Any, with $r = m$	Full row rank	Consistent for every $\mathbf{b}$
Any, with $r = n$	Full column rank	At most one solution for every $\mathbf{b}$

Numerical rank. In floating point, a nearly redundant equation leaves a tiny nonzero entry instead of an exact zero. Software therefore computes rank with a tolerance, usually from the singular value decomposition (Chapter 11). A matrix can be "rank 2" or "rank 1" depending on that tolerance, which is a warning sign of ill-conditioning.

import numpy as np

def classify(A, b):
    r = np.linalg.matrix_rank(A)
    r_aug = np.linalg.matrix_rank(np.column_stack([A, b]))
    n = A.shape[1]
    if r_aug > r:
        return "no solution"
    if r == n:
        return "exactly one solution"
    return f"infinitely many ({n - r} free parameters)"

A = np.array([[1., 2, -1, 3], [2, 4, 0, 2], [1, 2, 1, -1]])
print(classify(A, np.array([4., 6, 2])))        # infinitely many (2 free parameters)
print(classify(A, np.array([4., 6, 3])))        # no solution
B = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
print(classify(B, np.array([8., -11, -3])))     # exactly one solution

# Rank in floating point depends on a tolerance
C = np.array([[1., 1.], [1., 1. + 1e-12]])
print(np.linalg.matrix_rank(C), np.linalg.matrix_rank(C, tol=1e-9))   # 2 1

Key takeaways

• The rank is the number of pivots: the number of independent equations, and the dimension of the column space.
• Compare $r = \operatorname{rank} A$ with $r' = \operatorname{rank},[,A \mid \mathbf{b},]$: if $r' > r$ there is no solution; otherwise there are $n - r$ free parameters.
• $\operatorname{rank} A + \dim \operatorname{Nul} A = n$ (rank-nullity).
• Full row rank means solvable for every $\mathbf{b}$, and full column rank means at most one solution; numerically, rank depends on a tolerance.

🧪 Practice

1. Find the rank of $\begin{bmatrix} 1 & 2 & 3 \ 2 & 4 & 6 \ 1 & 0 & 1 \end{bmatrix}$ and the dimension of its null space.
2. A $4 \times 6$ matrix $A$ has rank 4. How many free variables does $A\mathbf{x} = \mathbf{b}$ have, and is it consistent for every $\mathbf{b}$?
3. For which values of $k$ does $x + y + kz = 1$, $\ x + ky + z = 1$, $\ kx + y + z = 1$ have no solution, exactly one solution, or infinitely many?
4. Interview: Can a $5 \times 3$ matrix have rank 4? Can a $3 \times 5$ matrix $A$ make $A\mathbf{x} = \mathbf{b}$ uniquely solvable for some $\mathbf{b}$? Justify both answers in one sentence each. Hint: bound the number of pivots by the number of rows and by the number of columns.

<a id="35-applications-of-linear-systems"></a>

3.5 Applications of Linear Systems

Linear systems appear wherever quantities are conserved, balanced, or combined in fixed proportions; this subchapter models five such situations and solves each one with the tools of this chapter.

Network and Traffic Flow

Why it matters. Road grids, pipe systems, power grids, data networks, and payment flows all obey a conservation law: at every junction, what flows in must flow out. Each junction contributes one linear equation. Free variables represent genuine routing choices, and the solution set describes every feasible way to route the flow.

Model. A network has nodes (junctions) connected by directed edges (streets, pipes) carrying unknown flows $x_1, \dots, x_n$. Some edges carry known external flow into or out of the network. At every node:

$$\text{total flow in} = \text{total flow out}.$$

For the whole network, total external inflow must equal total external outflow.

Example. Four one-way streets form a loop around a city block. External traffic (vehicles per hour) enters and leaves at the corners:

                 400
                  |
                  v
              +-------+        x1         +-------+
              |   A   |------------------>|   B   |------> 300
              +-------+                   +-------+
                  ^                           |
                  | x4                        | x2
                  |                           v
              +-------+        x3         +-------+
     300 <----|   D   |<------------------|   C   |<------ 200
              +-------+                   +-------+

Node	Flow in	Flow out	Equation
A	$400 + x_4$	$x_1$	$x_1 - x_4 = 400$
B	$x_1$	$x_2 + 300$	$x_1 - x_2 = 300$
C	$x_2 + 200$	$x_3$	$x_2 - x_3 = -200$
D	$x_3$	$x_4 + 300$	$x_3 - x_4 = 300$

Total inflow is $400 + 200 = 600$ and total outflow is $300 + 300 = 600$, so the network is balanced. Row reduction gives one free variable:

$$\mathbf{x} = (400, 100, 300, 0) + t,(1, 1, 1, 1), \qquad t = x_4.$$

The parameter $t$ is traffic that circulates around the block: it adds the same amount to every street without changing any node's balance. The direction $(1, 1, 1, 1)$ spans the null space of the coefficient matrix.

Redundancy. Adding all four node equations gives $0 = 0$, which is just the global balance. So one node equation is always redundant, and the rank is 3, not 4.

Physical constraints. Real flows must be nonnegative and below street capacity. Here $x_4 = t \ge 0$ forces $x_2 = 100 + t \ge 100$. Street B to C can never carry fewer than 100 vehicles per hour. Feasibility questions like this lead to linear programming.

import sympy as sp

x1, x2, x3, x4 = sp.symbols("x1:5")
node_eqs = [
    sp.Eq(400 + x4, x1),          # A: in = out
    sp.Eq(x1, x2 + 300),          # B
    sp.Eq(x2 + 200, x3),          # C
    sp.Eq(x3, x4 + 300),          # D
]
(sol,) = sp.linsolve(node_eqs, [x1, x2, x3, x4])
print(sol)                        # (x4 + 400, x4 + 100, x4 + 300, x4): one free variable

A, b = sp.linear_eq_to_matrix(node_eqs, [x1, x2, x3, x4])
print(A.rank())                   # 3: one node equation is redundant
print(A.nullspace())              # [Matrix([[1], [1], [1], [1]])]: circulation around the block

# Flows must be nonnegative, so x4 >= 0; the smallest feasible flow on each street is at x4 = 0
print([e.subs(x4, 0) for e in sol])   # [400, 100, 300, 0]

Key takeaways

• Each node contributes the linear equation "flow in = flow out".
• In a connected network, one node equation is redundant because the equations sum to the global balance.
• Free variables are routing choices; null-space directions such as $(1, 1, 1, 1)$ are circulations around loops.
• Physical limits (nonnegativity, capacity) restrict the parameters further.

🧪 Practice

1. In the example, the city closes street D to A ($x_4 = 0$). Find the flows on the other streets.
2. Suppose the 200 entering at C rises to 250 and the 300 leaving at D rises to 350. Re-solve the network and interpret the change.
3. Each street in the example has capacity 500 vehicles per hour. What range of values of $t$ is feasible?
4. Interview: For a connected network with $N$ nodes and $E$ edges, how many free parameters does the flow-conservation system have when it is consistent? Hint: the rank of the node equations is $N - 1$; the free parameters count independent loops.

Electrical Circuits and Kirchhoff's Laws

Why it matters. Circuit simulators such as SPICE turn a circuit into a large linear system and solve it, often millions of times during a simulation. For resistor networks with constant voltage sources, every governing law is linear, so circuit analysis is a direct application of this chapter.

The laws.

Law	Statement	Equation type
Ohm's law	The voltage drop across a resistor is $V = IR$	Relates current to voltage
Kirchhoff's current law (KCL)	At each junction, total current in = total current out	Conservation of charge, like network flow
Kirchhoff's voltage law (KVL)	Around any closed loop, the sum of voltage rises equals the sum of voltage drops	Conservation of energy

Sign convention. Choose a direction for each unknown current before writing equations. If a computed current is negative, it flows opposite to the chosen direction. No guesswork is needed.

Example. Two batteries feed a three-resistor network (resistances in ohms, voltages in volts):

               I1 -->                       I3 -->
      +-------[ R1 = 4 ]-------+-------[ R3 = 2 ]-------+
      |                        | P                      |
      |                        |                        |
  [ E1 = 20 V ]           [ R2 = 3 ]  | I2         [ E2 = 10 V ]
  (+ terminal up)              |      v            (+ terminal up)
      |                        |                        |
      +------------------------+------------------------+
                               Q

• KCL at P: $I_1 = I_2 + I_3$. The equation at Q is the same, so it is redundant.
• KVL, left loop (up through $E_1$, across $R_1$, down through $R_2$): $20 = 4I_1 + 3I_2$.
• KVL, right loop (up through $R_2$ against $I_2$, across $R_3$ with $I_3$, down through $E_2$ from $+$ to $-$): $3I_2 - 2I_3 - 10 = 0$.

$$ \left[\begin{array}{ccc|c} 1 & -1 & -1 & 0 \ 4 & 3 & 0 & 20 \ 0 & 3 & -2 & 10 \end{array}\right] \ \Longrightarrow
I_1 = \tfrac{35}{13} \approx 2.69, \quad I_2 = \tfrac{40}{13} \approx 3.08, \quad I_3 = -\tfrac{5}{13} \approx -0.38. $$

The negative $I_3$ means that current in $R_3$ flows from right to left: the 10 V source also pushes current toward P. A quick check on the left loop: $4 \cdot \tfrac{35}{13} + 3 \cdot \tfrac{40}{13} = \tfrac{260}{13} = 20$.

At scale. Nodal analysis chooses node voltages as unknowns and writes KCL at every node but one. For resistor networks this produces a symmetric coefficient matrix $G$, the conductance matrix, and $G\mathbf{v} = \mathbf{i}$ has a unique solution once one node is grounded. Chapter 10 explains why such matrices behave so well.

import sympy as sp

I1, I2, I3 = sp.symbols("I1 I2 I3")
eqs = [
    sp.Eq(I1, I2 + I3),             # KCL at node P
    sp.Eq(4*I1 + 3*I2, 20),         # KVL, left loop: source voltage = sum of IR drops
    sp.Eq(3*I2 - 2*I3, 10),         # KVL, right loop, with E2's + terminal up
]
sol = sp.solve(eqs, [I1, I2, I3])
print(sol)                                         # {I1: 35/13, I2: 40/13, I3: -5/13}
print({k: round(float(v), 3) for k, v in sol.items()})

# Power check: power delivered by the sources equals power dissipated in the resistors
delivered = 20 * sol[I1] + 10 * (-sol[I3])          # E2 drives current -I3 upward through itself
dissipated = 4 * sol[I1]**2 + 3 * sol[I2]**2 + 2 * sol[I3]**2
print(delivered == dissipated)                      # True

Key takeaways

• Ohm's law and Kirchhoff's current and voltage laws are all linear in the unknown currents.
• KCL gives one equation per junction, with one of them redundant; KVL gives one equation per independent loop.
• Choose current directions first; negative results indicate flow in the opposite direction.
• Large circuits are solved through nodal analysis, which leads to symmetric linear systems.

🧪 Practice

1. Re-solve the example with $E_2$ reversed (its $+$ terminal down). Which currents change sign or size?
2. Replace $R_3$ by a wire ($R_3 = 0$). Set up and solve the new system. What happens to $I_2$, and why?
3. Three resistors $R_1, R_2, R_3$ are connected in parallel across a single source $E$. Write KCL and KVL equations and show that the total current is $E,(1/R_1 + 1/R_2 + 1/R_3)$.
4. Interview: A circuit with $N$ nodes produces a nodal system whose coefficient matrix is singular until you ground one node. Why is the matrix singular, and why does grounding fix it? Hint: adding the same constant to every node voltage changes no current, so what vector lies in the null space?

Chemical Equation Balancing

Why it matters. Atoms are neither created nor destroyed in a chemical reaction, so each element must appear in equal numbers on both sides. The unknown coefficients satisfy a homogeneous linear system, and its null space gives the balanced equation. The method works automatically for reactions too tangled to balance by inspection.

Model. Represent each molecule by its composition vector, the count of each element in it. For propane combustion,

$$x_1,\mathrm{C_3H_8} + x_2,\mathrm{O_2} ;\to; x_3,\mathrm{CO_2} + x_4,\mathrm{H_2O},$$

use the element order (C, H, O):

Element	$\mathrm{C_3H_8}$	$\mathrm{O_2}$	$\mathrm{CO_2}$	$\mathrm{H_2O}$
C	3	0	1	0
H	8	0	0	2
O	0	2	2	1

Conservation requires $x_1\mathbf{r}_1 + x_2\mathbf{r}_2 = x_3\mathbf{p}_1 + x_4\mathbf{p}_2$. Moving the products to the left gives the homogeneous system

$$ \begin{bmatrix} 3 & 0 & -1 & 0 \ 8 & 0 & 0 & -2 \ 0 & 2 & -2 & -1 \end{bmatrix} \begin{bmatrix} x_1 \ x_2 \ x_3 \ x_4 \end{bmatrix} = \mathbf{0}. $$

Reduction leaves $x_4$ free, with $x_1 = \tfrac14 x_4$, $x_2 = \tfrac54 x_4$, $x_3 = \tfrac34 x_4$. The smallest positive integer solution, at $x_4 = 4$, is

$$\mathrm{C_3H_8} + 5,\mathrm{O_2} ;\to; 3,\mathrm{CO_2} + 4,\mathrm{H_2O}.$$

Reading the null space.

$\dim \operatorname{Nul} A$	Meaning
0	Only the trivial solution: the reaction as written is impossible (wrong formula or missing species)
1	Unique up to scaling; scale to the smallest positive integers
2 or more	Several independent reactions are mixed together; chemistry, not algebra, decides the ratios

Positivity. A null-space vector can have mixed signs. A valid balancing needs all coefficients positive. If no positive vector exists in the null space, some species is on the wrong side of the arrow.

import re
import sympy as sp

def composition(formula):
    """Element counts for simple formulas without parentheses, e.g. 'C3H8'."""
    counts = {}
    for elem, num in re.findall(r"([A-Z][a-z]?)(\d*)", formula):
        counts[elem] = counts.get(elem, 0) + int(num or 1)
    return counts

def balance(reactants, products):
    molecules = reactants + products
    comps = [composition(m) for m in molecules]
    elements = sorted(set().union(*comps))
    sign = [1] * len(reactants) + [-1] * len(products)      # products move left with a minus sign
    A = sp.Matrix([[s * c.get(e, 0) for c, s in zip(comps, sign)] for e in elements])
    basis = A.nullspace()
    if len(basis) != 1:
        raise ValueError(f"null space has dimension {len(basis)}: no unique balancing")
    v = basis[0] * sp.ilcm(*[entry.q for entry in basis[0]])  # clear denominators
    v = v / sp.igcd(*v)                                       # reduce to the smallest integers
    if any(entry < 0 for entry in v):
        v = -v
    if any(entry <= 0 for entry in v):
        raise ValueError("no balancing with all coefficients positive")
    return dict(zip(molecules, v))

print(balance(["C3H8", "O2"], ["CO2", "H2O"]))
# {'C3H8': 1, 'O2': 5, 'CO2': 3, 'H2O': 4}
print(balance(["KMnO4", "HCl"], ["KCl", "MnCl2", "Cl2", "H2O"]))
# {'KMnO4': 2, 'HCl': 16, 'KCl': 2, 'MnCl2': 2, 'Cl2': 5, 'H2O': 8}

Key takeaways

• Conservation of each element gives one homogeneous linear equation in the unknown coefficients.
• Columns are composition vectors, with products entering with a minus sign.
• A one-dimensional null space gives a unique balancing up to scale; clear denominators to get the smallest integers.
• A trivial null space, or one with no positive vector, signals an impossible reaction; a larger null space signals several independent reactions.

🧪 Practice

1. Balance $\mathrm{Al} + \mathrm{O_2} \to \mathrm{Al_2O_3}$ by setting up and solving the homogeneous system.
2. Balance $\mathrm{C_6H_{12}O_6} + \mathrm{O_2} \to \mathrm{CO_2} + \mathrm{H_2O}$ (cellular respiration).
3. Try to balance $\mathrm{H_2} + \mathrm{O_2} \to \mathrm{H_2O} + \mathrm{H_2O_2}$. What is the dimension of the null space, and what does it mean chemically?
4. Interview: Your balancing code returns a null-space vector with mixed signs, such as $(2, -1, 3)$. What does that tell you about the reaction as written? Hint: a negative coefficient on a reactant means that species is actually produced.

Polynomial Interpolation

Why it matters. Finding the polynomial that passes exactly through given data points is a building block of numerical integration, ODE solvers, computer graphics, and Shamir's secret-sharing scheme, which works over the finite fields of Section 1.2. The problem is linear even though polynomials are not: the unknowns are the coefficients, and every data point gives one linear equation in them.

Setup. Given $n + 1$ points $(t_0, y_0), \dots, (t_n, y_n)$ with distinct $t_i$, find $p(t) = c_0 + c_1t + \cdots + c_nt^n$ with $p(t_i) = y_i$. Each condition reads

$$c_0 + c_1 t_i + c_2 t_i^2 + \cdots + c_n t_i^n = y_i,$$

which is linear in $c_0, \dots, c_n$. Together they form $V\mathbf{c} = \mathbf{y}$ with the Vandermonde matrix

$$V = \begin{bmatrix} 1 & t_0 & t_0^2 & \cdots & t_0^n \ 1 & t_1 & t_1^2 & \cdots & t_1^n \ \vdots & & & & \vdots \ 1 & t_n & t_n^2 & \cdots & t_n^n \end{bmatrix}.$$

Theorem. If the nodes $t_0, \dots, t_n$ are distinct, there is exactly one interpolating polynomial of degree at most $n$.

Proof. $V$ is square. If $V\mathbf{c} = \mathbf{0}$, the polynomial with coefficients $\mathbf{c}$ has degree at most $n$ and $n + 1$ distinct roots, so it is the zero polynomial and $\mathbf{c} = \mathbf{0}$. Then $\operatorname{Nul} V = {\mathbf{0}}$, so every column is a pivot column, the rank is $n + 1$, and $V\mathbf{c} = \mathbf{y}$ has exactly one solution for every $\mathbf{y}$ (Section 3.4). $\blacksquare$

Example. Through $(1, 2)$, $(2, 3)$, $(3, 6)$:

$$ \left[\begin{array}{ccc|c} 1 & 1 & 1 & 2 \ 1 & 2 & 4 & 3 \ 1 & 3 & 9 & 6 \end{array}\right] \ \xrightarrow{\ \text{rref}\ }
\left[\begin{array}{ccc|c} 1 & 0 & 0 & 3 \ 0 & 1 & 0 & -2 \ 0 & 0 & 1 & 1 \end{array}\right] \ \Rightarrow\ p(t) = 3 - 2t + t^2. $$

   y
   6 |                 *  (3, 6)
   5 |               .
   4 |              .
   3 |           *  (2, 3)               p(t) = 3 - 2t + t^2
   2 |     *  .  (1, 2)                  passes through all three points
   1 |
     +-----+-----+-----+--- t
           1     2     3

A numerical warning. For many equally spaced nodes, the Vandermonde matrix is severely ill-conditioned, and the interpolant itself can oscillate wildly near the ends of the interval (Runge's phenomenon). In practice, use Chebyshev nodes and the Lagrange or Newton form of the polynomial, or replace one high-degree polynomial with piecewise low-degree splines.

import numpy as np

t = np.array([1., 2, 3])
y = np.array([2., 3, 6])
V = np.vander(t, increasing=True)          # columns: 1, t, t^2
c = np.linalg.solve(V, y)
print(c)                                   # approximately [ 3. -2.  1.]  ->  p(t) = 3 - 2t + t^2
print(np.polyval(c[::-1], 4.0))            # about 11.0: the interpolant evaluated at a new point

# Equally spaced nodes make V badly conditioned as the degree grows
for n in (5, 10, 20):
    nodes = np.linspace(0, 1, n + 1)
    print(n, f"{np.linalg.cond(np.vander(nodes, increasing=True)):.1e}")

Key takeaways

• Interpolation conditions $p(t_i) = y_i$ are linear in the polynomial's coefficients, giving a Vandermonde system $V\mathbf{c} = \mathbf{y}$.
• With $n + 1$ distinct nodes, there is exactly one interpolating polynomial of degree at most $n$.
• The proof is a null-space argument: a nonzero polynomial of degree at most $n$ cannot have $n + 1$ roots.
• High-degree interpolation on equally spaced nodes is numerically fragile; prefer better nodes, better bases, or splines.

🧪 Practice

1. Find the polynomial of degree at most 2 through $(0, 1)$, $(1, 0)$, and $(2, 3)$.
2. Find the cubic through $(-1, 0)$, $(0, 1)$, $(1, 0)$, $(2, 3)$. Is it actually of degree 3?
3. Explain what goes wrong with the Vandermonde system if two nodes coincide, $t_i = t_j$ with $i \ne j$, in terms of rank and consistency.
4. Interview: In Shamir's secret sharing, a secret $s$ is stored as $p(0)$ for a random polynomial $p$ of degree $k - 1$ over a finite field, and each participant receives one point $(t_i, p(t_i))$. Why do any $k$ participants recover the secret, while $k - 1$ learn nothing about it? Hint: count equations and unknowns, and for $k - 1$ points, describe the solution set as $s$ varies.

Leontief Input-Output Models

Why it matters. Wassily Leontief won the 1973 Nobel Prize in economics for modeling an economy as sectors that consume each other's output. The central question is how much each sector must produce to meet outside demand when production itself uses up inputs. The answer comes from a single linear system, and the model is still used to plan production and trace supply-chain effects.

Model. Suppose the economy has $n$ sectors.

• $x_j$ is the total output of sector $j$, in dollars.
• The consumption matrix $C$ has entry $c_{ij}$ equal to the dollar amount of sector $i$'s output needed to produce one dollar of sector $j$'s output. Column $j$ is sector $j$'s "recipe" of inputs.
• $C\mathbf{x}$ is the intermediate demand, the output consumed by the sectors themselves.
• $\mathbf{d}$ is the final demand from households, government, and exports.

Production must cover both:

$$\mathbf{x} = C\mathbf{x} + \mathbf{d} \quad\Longleftrightarrow\quad (I - C),\mathbf{x} = \mathbf{d},$$

where $I$ is the identity matrix (ones on the diagonal, so that $I\mathbf{x} = \mathbf{x}$; Section 4.2).

   total output  x   =     C x        +        d
                       (used up by the      (delivered to households,
                        sectors to make      government, exports)
                        their own output)

Example. Three sectors (manufacturing, agriculture, services), with

$$C = \begin{bmatrix} 0.50 & 0.40 & 0.20 \ 0.20 & 0.30 & 0.10 \ 0.10 & 0.10 & 0.30 \end{bmatrix}, \qquad \mathbf{d} = \begin{bmatrix} 50 \ 30 \ 20 \end{bmatrix}.$$

Column 1 says that one dollar of manufacturing output consumes 50 cents of manufacturing, 20 cents of agriculture, and 10 cents of services. Solving $(I - C)\mathbf{x} = \mathbf{d}$ gives $\mathbf{x} \approx (225.9,\ 118.5,\ 77.8)$. Meeting a final demand of 100 requires a total output of about 422, because of the internal demand.

Why a sensible answer exists. If every column of $C$ sums to less than 1, each sector needs less than one dollar of inputs per dollar of output, and the economy is productive. Then $(I - C)\mathbf{x} = \mathbf{d}$ has exactly one solution for every $\mathbf{d}$, and that solution has no negative entries whenever $\mathbf{d}$ has none. The reason is the expansion

$$\mathbf{x} = \mathbf{d} + C\mathbf{d} + C^2\mathbf{d} + C^3\mathbf{d} + \cdots$$

Read it as rounds of production. First produce $\mathbf{d}$. That needs inputs $C\mathbf{d}$, which need inputs $C^2\mathbf{d}$, and so on. Each round is smaller than the last, so the sum converges. Section 13.2 makes the conditions precise with Perron-Frobenius theory.

import numpy as np

C = np.array([[0.50, 0.40, 0.20],      # column j: inputs consumed per dollar of output of sector j
              [0.20, 0.30, 0.10],      # rows and columns: manufacturing, agriculture, services
              [0.10, 0.10, 0.30]])
d = np.array([50., 30, 20])            # final demand

x = np.linalg.solve(np.eye(3) - C, d)  # (I - C) x = d
print(x.round(1))                      # [225.9 118.5  77.8]
print(np.allclose(x, C @ x + d))       # True: output covers internal use plus final demand

# Rounds of production: d, then the inputs C d needed to make it, then C(C d), ...
total, round_ = np.zeros(3), d.copy()
for _ in range(200):
    total += round_
    round_ = C @ round_
print(total.round(1))                  # [225.9 118.5  77.8]: the same answer

Key takeaways

• In the Leontief model, column $j$ of $C$ lists the inputs needed per dollar of sector $j$'s output.
• Output must cover intermediate demand plus final demand: $\mathbf{x} = C\mathbf{x} + \mathbf{d}$, that is, $(I - C)\mathbf{x} = \mathbf{d}$.
• When every column sum of $C$ is below 1, the system has a unique nonnegative solution for every nonnegative $\mathbf{d}$.
• The solution expands as rounds of production $\mathbf{d} + C\mathbf{d} + C^2\mathbf{d} + \cdots$, which explains why total output exceeds final demand.

🧪 Practice

1. For a two-sector economy with $C = \begin{bmatrix} 0.2 & 0.3 \ 0.4 & 0.1 \end{bmatrix}$ and $\mathbf{d} = (100, 50)$, solve $(I - C)\mathbf{x} = \mathbf{d}$.
2. In the three-sector example, final demand for services rises by 10. Compute the new production vector and the change in each sector's output.
3. Compute the first three rounds $\mathbf{d}$, $C\mathbf{d}$, $C^2\mathbf{d}$ for the two-sector economy in Exercise 1 and compare their running sum with the exact solution.
4. Interview: An analyst's consumption matrix has one column summing to 1.2. What does that mean economically, and what can go wrong when solving $(I - C)\mathbf{x} = \mathbf{d}$? Hint: that sector consumes more than a dollar of inputs per dollar of output; think about what happens to the rounds $C^k\mathbf{d}$.

<a id="4-matrix-algebra"></a>