The word fangcheng, the modern Chinese term for “equation,” is older, and carries more meaning, than most people realize.
Chapter 8 of The Nine Chapters on the Mathematical Art is titled Fangcheng (“Rectangular Arrays”); the book was compiled and edited over a long period of accumulation. Its method lays out counting rods as an array of numbers: the coefficients of each equation are arranged along a vertical column, the unknowns are eliminated step by step by multiplying and subtracting, and the solution is then found by back substitution. This is an important record of the ancient idea of elimination; what was used at the time was counting rods on a counting board, not the modern abacus.
In modern notation, writing each equation horizontally as a row gives the augmented matrix, and the corresponding elimination steps are written as row operations. The ancient arrays and modern matrices follow different conventions of orientation, but both use equivalence transformations that leave the solutions unchanged. This connection helps us understand the history, but it does not mean that the ancients already used the modern concept of a matrix.
By the middle of the nineteenth century, systems of linear equations had begun to play an increasingly important role in problems of physics and engineering.
In 1845, Gustav Kirchhoff studied the relations between currents and voltages in electrical networks. Conservation of current at the nodes and the voltage relations around loops, together with the constitutive equations of the components, provided a systematic framework for circuit analysis. For a linear resistive network these relations form a system of linear equations; the number of independent equations must take the topological constraints of the network into account, and is not obtained by simply adding up the numbers of nodes, branches, or loops.
Under a suitable linear model, the connection and equilibrium conditions of a physical system can be translated into a system of linear equations. Nonlinear components, by contrast, usually require nonlinear solves and local linearization.
This framework extends to modern circuit simulation: many algorithms solve linear systems over and over in the course of their iterations. Multi-antenna channel equalization in communications also often involves linear equations and least-squares problems. Different applications rest on different modeling assumptions, yet they share one body of linear algebra tools.
More profoundly, solving a system of equations means far more than arithmetic elimination. From the modern geometric viewpoint, solving Ax=b is at once a search for the common intersection point of several hyperplanes in a high-dimensional space and a test of whether the target vector b can be formed as a linear combination of the column vectors of the matrix.
The grain calculations on ancient counting boards, the resistive networks in the notes of nineteenth-century physicists, and today’s simulations of electromagnetic fields on supercomputers can all, under suitable models, be reduced to systems of linear equations. This consistency across civilizations, engineering, and pure geometry is an important reason why systems of linear equations sit at the heart of linear algebra.
The elimination algorithm of The Nine Chapters is intuitive and effective, but what structure lies behind it? Why can those operations be applied freely without changing the solutions of the system? §6.1 starts from this question and interprets row operations as invertible linear transformations—left multiplication by elementary matrices—thereby raising the arithmetic act of “solving equations” to the language of linear algebra.
Each of these elementary operations corresponds to an invertible elementary matrix, so they preserve the dimensions of the row space and the column space, which guarantees the legitimacy of every subsequent reduction step in the chapter.
§6.2 introduces the row echelon form and the reduced row echelon form, the end point of elimination: once a matrix has been brought to this canonical form, the structure of the solutions—no solution, a unique solution, or infinitely many solutions—becomes clear at a glance. Here the concept of rank receives its most concrete computational definition, and it formally connects with the concept of dimension from Chapter 4 on abstract linear spaces.
§6.3 contains the central algorithms of the chapter: Gaussian elimination and the LU decomposition. The former is the procedure for solving; the latter is the same procedure in the language of matrices—factoring A as the product of a lower triangular matrix L and an upper triangular matrix U. The significance of the LU decomposition goes beyond computational efficiency: it records every step of the elimination permanently in L, so that the same coefficient matrix can be reused for different right-hand-side vectors, which is one of the basic design principles of modern numerical linear algebra.
§6.4 closes with three concrete applications: the balance equations of network flows, the conservation of atoms in chemical reactions, and the static equilibrium of mechanical systems. These three problems come from different fields, yet mathematically they are all Ax=b—they share one theory and one set of algorithms.
Once you have read this chapter, you will see a different depth of field in the ancient name of the fangcheng method.
We begin with the most basic matrix operations and set out to explore the mysteries of systems of linear equations.
Elementary matrix operations are the basic tools for solving systems of linear equations. In essence they are an algebraic upgrade of the three operations for solving equations that you already know from high school—multiplying by a constant, interchanging two equations, and adding two equations. Let us first look at a concrete example to get a feel for this upgrade:
Multiplying row 0 of the system by 2 is completely equivalent to “left-multiplying by the diagonal matrix E0(2)”:
E0(2)[2001]A[1235]=[2265]=result of row scaling[2⋅row0(A)row1(A)]
This observation reveals the central theme of this section: every elimination operation done by hand corresponds to left multiplication by a particular elementary matrix. Once you master this correspondence, Gaussian elimination turns from a “technique” into a “structure.”
Row scaling means multiplying one row of a matrix by a nonzero constant.
Row scaling does not change the dimension of the row space, because it merely stretches one row vector by a nonzero factor and does not change linear dependence. At the same time, the dimension of the column space is also unchanged, because Ei(c) is invertible and an invertible transformation does not change the dimension of the image. Since the rank is defined as the dimension of the column space, row scaling does not change the rank of the matrix.
In a system of linear equations Ax=b, row scaling amounts to multiplying both sides of one equation by a nonzero constant, which does not change the solution set of the equation.
Row exchange means interchanging the positions of two rows of a matrix.
Row exchange does not change the dimension of the row space, because interchanging two rows merely reorders the row vectors without changing what their linear combinations can produce. Since Pij is an invertible matrix (in fact an orthogonal matrix), the dimension of the column space is also unchanged. Therefore row exchange preserves the rank of the matrix.
The row exchange operation amounts to reordering the equations of a system of linear equations, which does not change the solution set of the system.
Row elimination means adding a multiple of one row to another row, with the aim of eliminating the entry in a particular position.
Row elimination does not change the dimension of the row space, because it merely replaces one row vector by a linear combination of that vector and another row vector, and the span of all the row vectors stays the same. Since Eij(c) is an invertible matrix (its inverse is Eij(−c)), the dimension of the column space is also unchanged. Therefore the row elimination operation preserves the rank of the matrix.
Row elimination is the central operation for solving systems of linear equations: it simplifies the system by eliminating variables, but it does not change the solution set of the system.
6.1.4 The Relation Between Row Operations and Column Operations¶
Column operations and row operations are dual to each other, and the connection between them is made through the matrix transpose. A column operation on the matrix A is equivalent to the corresponding row operation on its transpose A⊤.
From the property above and the fact that row operations preserve the dimension of the column space, we can deduce that column operations preserve the dimension of the row space. Since dim(Col(A))=dim(Row(A⊤)), column operations do not change the dimension of the row space, and row operations do not change the dimension of the column space; hence all elementary operations preserve the rank of the matrix.
Column scaling
Column scaling means multiplying one column of a matrix by a nonzero constant; it is equivalent to right multiplication by a diagonal elementary matrix Ej(c), written AEj(c).
Column scaling does not change the dimension of the column space, because it merely stretches one column vector by a nonzero factor and does not change linear dependence. By the definition of the rank of a matrix (the dimension of the column space), column scaling does not change the rank of the matrix.
Column exchange
Column exchange means interchanging the positions of two columns of a matrix; it is equivalent to right multiplication by a permutation matrix Pij, written APij.
Column exchange does not change the dimension of the column space, because interchanging two columns merely reorders the column vectors without changing what their linear combinations can produce. Therefore column exchange preserves the rank of the matrix.
Column elimination
Column elimination means adding a multiple of one column to another column; it is equivalent to right multiplication by an elementary matrix Eji(c), written AEji(c).
Column elimination does not change the dimension of the column space, because it merely replaces one column vector by a linear combination of that vector and another column vector, and the span of all the column vectors stays the same. Therefore the column elimination operation preserves the rank of the matrix.
Every elementary matrix operation preserves the dimensions of the row space and the column space of the matrix; this is the key to understanding the structure of the solutions of a system of linear equations. We give a slightly more general statement.
The core property of elementary matrix operations is that they preserve the dimensions of the row space and the column space, and hence the rank of the matrix. This property is the foundation of Gaussian elimination and the other algorithms for solving systems of linear equations. Row operations apply an invertible linear transformation on the output side and leave the solution set unchanged, which makes them the main tool for solving systems of linear equations; column operations have important applications in matrix decompositions and in understanding the structure of the solution space. Understanding how elementary operations affect the subspace structure of a matrix is the key to mastering the core ideas of linear algebra.
6.2 Reduced Row Echelon Form and the Geometric Structure of Solutions¶
The reduced row echelon form is not only a key intermediate step in solving systems of linear equations; it also gives an intuitive expression of the geometric structure of the solutions. Through the reduced row echelon form we can clearly see the dimension and shape of the solution space of a system and its relation to the subspaces of the coefficient matrix.
6.2.1 Square Upper Triangular Matrices with Nonzero Diagonal and Their Geometric Interpretation¶
An upper triangular matrix with nonzero diagonal entries is like a “staircase” system of equations: each equation is “simpler” than the one before it, so we can start from the last equation and solve backward step by step. In §5.3.3, Theorem 3, we already proved that such a system necessarily has a solution; now we work through an actual example and solve it.
Why do upper triangular matrices have this special structure?
From the viewpoint of transformations, an upper triangular matrix has a kind of “hierarchy”:
It maps the standard basis vectors {e0,e1,…,en−1} to an image with a special structure
This image can be decomposed into nested, properly contained subspaces: {0}⊂span{e0}⊂span{e0,e1}⊂⋯⊂Rn. This is because
every diagonal entry is nonzero, which guarantees that every standard basis vector makes a “new contribution” to the image
the matrix has full rank (rank = n), so the image is the whole of Rn
each additional basis vector increases the dimension of the image by 1
An intuitive analogy. It is like the floors of a building: the upper floors can affect the lower floors, but the lower floors cannot affect the upper ones. This “one-way dependence” lets us start from the “top floor” and solve downward floor by floor.
This structure is especially important in numerical computation, because it guarantees the stability and efficiency of the solution process.
6.2.2 Row Echelon Matrices and the Reduced Row Echelon Form¶
A square upper triangular matrix whose diagonal entries are all nonzero is a REF whose pivots all lie on the diagonal; a general upper triangular matrix need not be a REF.
A matrix in row echelon form (Row Echelon Form, REF) is a special kind of upper triangular matrix (conversely, an upper triangular matrix need not be in row echelon form); its image can be decomposed into nested subspaces, but the containments are not necessarily proper.
Let us next consider a simpler form, which can be converted to and from a row echelon matrix by row scaling and row elimination: the reduced row echelon matrix (Reduced Row Echelon Form, RREF), which simplifies this geometric representation further.
6.2.3 Reading the Rank and the Subspace Dimensions from the Block Structure of the RREF¶
§6.2.1 showed that an upper triangular matrix with nonzero diagonal has a nested “flag” structure of subspaces: each new pivot contributes an entirely new direction to the image, and back substitution solves step by step from the last row upward. The RREF and the upper triangular matrix share essentially the same “staircase of pivots” geometry; the only difference is that an upper triangular matrix allows nonzero entries above the pivots, while the RREF further requires the entries above and below each pivot to be cleared to zero and each pivot to equal 1. This extra reduction costs very little (row elimination plus row scaling), but it brings a decisive benefit: the block structure of the RREF can be read off directly, with no back substitution needed.
However, the pivots of an RREF do not fall neatly on the diagonal as they do for a square upper triangular matrix. Take the following matrix as an example:
The pivots are in columns 0, 1, and 3, while columns 2 and 4 have no pivot. This phenomenon of “pivots skipping columns” tells us that the order of the variables needs to be adjusted before the RREF can be rearranged into a standard block form. Below we illustrate the whole process with a complete numerical example, and then state the theorem.
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ──────────────────────────────────────────────
def print_header(title):
print("=" * 60)
print(f" {title}")
print("=" * 60)
def print_step(step, desc):
print(f"\n▶ Step {step}: {desc}")
print("-" * 40)
def compare_print(label, actual, expected_desc):
print(f"\n[{label}]")
print(f" Computed value: {actual}")
print(f" Expected value: {expected_desc}")
# ============================================================
print_header("Example 6.1 | Solving a Homogeneous RREF System and Verifying the Kernel")
# ============================================================
# --- Step 1: define the RREF matrix R ---
print_step(1, "Define the RREF matrix R")
R = np.array([
[1, 0, 2, 0, 5],
[0, 1, 3, 0, 4],
[0, 0, 0, 1, -2],
[0, 0, 0, 0, 0]
], dtype=float)
print("R =\n", R)
# --- Step 2: read a basis of the kernel directly from the block structure ---
print_step(2, "Basis of the kernel (read directly from the free-variable columns)")
# pivot columns: 0, 1, 3 ; free columns: 2, 4
# assign unit values to the free variables (x₂, x₄) in turn
n0 = np.array([-2, -3, 1, 0, 0], dtype=float) # x₂ = 1, x₄ = 0
n1 = np.array([-5, -4, 0, 2, 1], dtype=float) # x₂ = 0, x₄ = 1
print("n₀ =", n0)
print("n₁ =", n1)
# --- Step 3: check R @ nᵢ = 0 ---
print_step(3, "Check R·nᵢ = 0 (nᵢ belongs to Null(R))")
compare_print("R · n₀", R @ n0, "zero vector [0, 0, 0, 0]")
compare_print("R · n₁", R @ n1, "zero vector [0, 0, 0, 0]")
# --- Step 4: check the dimensions with the rank–nullity theorem ---
print_step(4, "Check with the rank–nullity theorem")
rank_R = np.linalg.matrix_rank(R)
nullity_R = R.shape[1] - rank_R
compare_print("rank(R)", rank_R, "3 (three pivot columns)")
compare_print("nullity(R)", nullity_R, "2 (two free columns)")
print(f"\n rank(R) + nullity(R) = {rank_R} + {nullity_R} = {rank_R + nullity_R}"
f" = n = {R.shape[1]} ✓")
With this concrete example as a foundation, we can now state the general theorem.
6.2.4 Nonhomogeneous Systems of Linear Equations: Particular Solution plus Homogeneous Solution¶
§6.2.3 settled the homogeneous system Rx=0: the structure of the kernel is completely determined by the block matrix [IrF], and the basis matrix is N=P[−FIn−r].
Now consider the general case Ax=b, b=0. Compared with the homogeneous case, a nonhomogeneous system faces a new problem: a solution need not exist. Before discussing how to solve, we must first ask: what condition must b satisfy for the system to have a solution?
The hyperplane view. Each row of the system Ax=b is the equation of a hyperplane in Rn. Having a solution is equivalent to all the hyperplanes having a common point of intersection. If the RREF of the augmented matrix [R∣c] contains a row of the form
then that row requires “0=k,” which cannot be satisfied by any x; this shows that the hyperplanes contradict one another and have no common point of intersection.
The image view. Regard Ax=b as asking: does b lie in the column space (the image) Col(A) of A?
Because Ax is a linear combination of the columns of A, having a solution is equivalent to b being expressible as a linear combination of these columns. In terms of rank: b∈Col(A) if and only if appending b to A does not increase the dimension of the column space, that is,
The two views are equivalent: an inconsistent row appearing in the RREF is equivalent to b not lying in the image of A, which is equivalent to the rank of the augmented matrix exceeding the rank of the coefficient matrix.
The Structure of the Solutions When a Solution Exists: Particular Solution plus Homogeneous Solution¶
Once we have confirmed that the system has a solution, the structure of the solutions is described by the following theorem:
Let R=EA; the right-hand side must undergo the same row operations, giving c=Eb. Let bp be the first r components of c; consistency requires the remaining m−r components all to be zero. After the columns are rearranged we have yp=bp−Fyf. In this block structure the most natural choice of particular solution is to set all the free variables to zero (yf=0); then yp=bp, and restoring the order of the columns gives a particular solution. The complete solution is then obtained by adding to this particular solution the kernel already found in §6.2.3.
The three forms of the solutions are now clear at a glance:
Case
Condition (RREF)
Solution set
Dimension
No solution
An inconsistent row [0⋯0∣k], k=0, appears
Empty set
—
Unique solution
No inconsistent row, and no free variables (r=n)
A single point {xp}
0
Infinitely many solutions
No inconsistent row, and there are free variables (r<n)
Gaussian elimination is the systematic method for transforming the augmented matrix of a system of linear equations into row echelon form (REF) or reduced row echelon form (RREF). It provides the actual computational steps for reaching the matrix forms discussed in §6.2.
At the heart of Gaussian elimination is a sequence of elementary row operations that transform the augmented matrix into an equivalent row echelon matrix. These operations are in essence isomorphic transformations of the image, which guarantee that the solution set of the system is the same before and after the transformation.
6.3.2 Simplified Gaussian Elimination and the LU Decomposition¶
Before giving the definition, let us first settle a question of notation: the elimination process produces a product of elementary matrices
En−1⋯E1E0A=U,
and there are two possible choices of notation:
Core idea. Ideal Gaussian elimination brings the original matrix A step by step to an upper triangular (upper trapezoidal) matrix U using pure “row elimination operations” alone (without interchanging the order of the rows); if the pivots fall in turn on the diagonal and are all nonzero, U is also in row echelon form; otherwise, in general, it is only guaranteed to be upper triangular and need not be a REF.
The following code implements the algorithm above; please read it alongside the steps.
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ──────────────────────────────────────────────
def print_header(title):
print("=" * 60)
print(f" {title}")
print("=" * 60)
def print_step(step, desc):
print(f"\n▶ Step {step}: {desc}")
print("-" * 40)
def compare_print(label, actual, expected_desc):
print(f"\n[{label}]")
print(f" Computed value: {actual}")
print(f" Expected value: {expected_desc}")
# ──────────────────────────────────────────────
def simplified_lu_decomposition(A, tol=1e-12):
"""
Simplified LU decomposition (no row exchanges)
Parameters
----------
A : input matrix, of shape m×n
tol : numerical tolerance for deciding whether a pivot is zero
Returns
-------
L : lower triangular matrix, of shape m×m
U : upper triangular matrix, of shape m×n
success : Boolean indicating whether the decomposition succeeded
"""
m, n = A.shape
L = np.eye(m) # initialize L as the m×m identity matrix
U = A.copy().astype(float) # initialize U as a floating-point copy of A
for k in range(min(m - 1, n)):
# --- Step 1: check the pivot U[k, k] ---
if abs(U[k, k]) >= tol:
pass # the pivot is usable; go straight to elimination
else:
# --- Step 2: the pivot is close to zero; check below it in this column ---
column_below = U[k + 1:m, k]
if np.any(np.abs(column_below) >= tol):
# there is a nonzero entry below, but row exchanges are not allowed
print(f" [!] The pivot in row {k}, column {k} is zero, and there is a nonzero entry below it")
print(f" The simplified LU decomposition does not support row exchanges; the decomposition fails")
return L, U, False
else:
# the whole column is close to zero; skip this column
print(f" [i] The column from row {k}, column {k} down is close to zero; skipping")
continue
# --- Step 3: eliminate in column k ---
for i in range(k + 1, m):
if abs(U[i, k]) >= tol:
multiplier = U[i, k] / U[k, k] # elimination multiplier l_{ik}
L[i, k] = multiplier # store in L
U[i, k:n] = U[i, k:n] - multiplier * U[k, k:n]
U[i, k] = 0.0 # set to zero explicitly to remove floating-point residue
return L, U, True
# ============================================================
print_header("Example 6.2 | Simplified LU Decomposition")
# ============================================================
# --- Step 1: define the matrix A ---
print_step(1, "Define the matrix A")
A = np.array([
[ 2.0, 1.0, 1.0],
[ 4.0, -6.0, 0.0],
[-2.0, 7.0, 2.0]
], dtype=float)
print("A =\n", A)
# --- Step 2: perform the simplified LU decomposition ---
print_step(2, "Perform the simplified LU decomposition (no row exchanges)")
L, U, success = simplified_lu_decomposition(A)
if not success:
print("\n The decomposition failed; use the PLU decomposition with row exchanges instead.")
else:
print(" Decomposition succeeded ✓")
# --- Step 3: print L and U ---
print_step(3, "Result of the decomposition: L and U")
print("L =\n", L)
print("\nU =\n", U)
# --- Step 4: check L @ U = A ---
print_step(4, "Check L·U = A")
LU = L @ U
error = np.max(np.abs(A - LU))
compare_print("L · U", LU, "should recover the original matrix A")
compare_print("Maximum reconstruction error ‖A − LU‖∞", f"{error:.2e}", "should be close to machine precision ~1e-16")
# --- Step 5: check the upper triangular structure of U (all zeros below the diagonal) ---
print_step(5, "Check that U is upper triangular: U[i, j] = 0 when i > j")
below_diag = np.max(np.abs(np.tril(U, k=-1)))
compare_print("Largest absolute value below the diagonal max|U[i, j]| (i > j)", f"{below_diag:.2e}", "should be 0")
assert error < 1e-10, "Reconstruction error too large: L·U does not recover A"
assert below_diag < 1e-12, "U is not upper triangular: nonzero entries remain below the diagonal"
# --- Step 6: test with a tall matrix (m > n): the last elimination position (column 1) must not be skipped ---
print_step(6, "Tall-matrix test: the 3×2 matrix A_tall")
A_tall = np.array([
[1.0, 0.0],
[0.0, 1.0],
[0.0, 2.0]
])
L_t, U_t, ok_t = simplified_lu_decomposition(A_tall)
print("U_tall =\n", U_t)
err_t = np.max(np.abs(A_tall - L_t @ U_t))
below_t = np.max(np.abs(np.tril(U_t, k=-1)))
compare_print("‖A_tall − LU‖∞ / largest value below the diagonal of U", f"{err_t:.2e} / {below_t:.2e}", "both should be 0")
assert ok_t and err_t < 1e-10 and below_t < 1e-12, "The LU decomposition of a tall matrix should both reconstruct A and be upper triangular"
Why Is the LU Decomposition Written A=LU Rather Than LA=U?¶
Each step of Gaussian elimination left-multipliesA by an elementary lower triangular matrixEk, clearing the entries below the pivot in column k:
So from the viewpoint of the elimination operations, the notation LA=U is not entirely wrong—it describes “applying a sequence of row operations to A to obtain U.”
But there are three fundamental reasons for writing it in the end as A=LU:
Reason 1: the inverses of the Ek have a more elegant structure
Each Ek is a unit lower triangular matrix, and its inverse Ek−1 is obtained simply by changing the signs of the elimination multipliers, with exactly the same structure. Let
Then L is precisely the lower triangular matrix obtained by entering all the elimination multipliers directly into their corresponding positions, with no extra computation. This is exactly what the code below demonstrates: multiplying the Ek−1in order of increasing index (E0−1E1−1⋯En−1−1, that is, in the reverse of the order in which the eliminations are performed) gives a product that is exactly the L assembled directly from the multipliers of each round; if they are multiplied in order of decreasing index (En−1−1⋯E0−1), the product is still a unit lower triangular matrix but is not equal to L.
Reason 2: A=LU is a decomposition of A, while LA=U is a transformation of A
The meaning of a decomposition is “to express A as the product of two matrices of simple structure,” for use in solving, analysis, and storage. The meaning of LA=U, by contrast, is “to apply operations to A”; it describes a process rather than a structure.
Reason 3: the applications are built around A=LU
When solving the linear system Ax=b, the decomposition gives
Both substitution steps use L and U directly rather than the sequence of Ek. The form A=LU makes this procedure natural and symmetric.
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ──────────────────────────────────────────────
def print_header(title):
print("=" * 60)
print(f" {title}")
print("=" * 60)
def print_step(step, desc):
print(f"\n▶ Step {step}: {desc}")
print("-" * 40)
def compare_print(label, actual, expected_desc):
print(f"\n[{label}]")
print(f" Computed value:\n{actual}")
print(f" Expected value: {expected_desc}")
# ============================================================
print_header("Example 6.3 | The Order of the Product of the Inverse Elementary Matrices E_k⁻¹ and L")
# ============================================================
n = 4 # size for the demonstration; adjust as you like
np.random.seed(42) # fix the random seed so the results are reproducible
# --- Step 1: generate n unit lower triangular matrices E_0⁻¹, …, E_{n-1}⁻¹ ---
# Note: the elimination matrix E_k in the text places the multipliers, "with their signs changed," below the pivot in column k;
# its inverse E_k⁻¹ places the elimination multipliers "unchanged (with positive sign)" in the same positions.
# Here we directly generate the E_k⁻¹ of each round (with positive multipliers), matching L = E_0⁻¹E_1⁻¹⋯E_{n-1}⁻¹.
print_step(1, f"Generate {n} unit lower triangular matrices E_0⁻¹, …, E_{{n-1}}⁻¹ (each holding the multipliers of round k)")
Einv_list = []
for k in range(n):
Einv = np.eye(n)
num_fill = n - k - 1 # number of multipliers to fill in below the pivot in column k
if num_fill > 0:
Einv[k+1:, k] = np.round(np.random.uniform(1, 10, num_fill), 2)
Einv_list.append(Einv)
print(f"\nE_{k}⁻¹ =\n", Einv)
# --- Step 2: multiply in order E_0⁻¹ @ E_1⁻¹ @ … @ E_{n-1}⁻¹ (increasing index) ---
print_step(2, "Multiply in order: E_0⁻¹ · E_1⁻¹ · … · E_{n-1}⁻¹ (increasing index)")
forward = Einv_list[0].copy()
for Einv in Einv_list[1:]:
forward = forward @ Einv
print("Product in order =\n", forward)
print("\n ➜ The product in order is a unit lower triangular matrix and is exactly the L of the LU decomposition (see the check in Step 4)")
# --- Step 3: multiply in reverse order E_{n-1}⁻¹ @ … @ E_1⁻¹ @ E_0⁻¹ (decreasing index) ---
print_step(3, "Multiply in reverse order: E_{n-1}⁻¹ · … · E_1⁻¹ · E_0⁻¹ (decreasing index)")
reverse = Einv_list[-1].copy()
for Einv in reversed(Einv_list[:-1]):
reverse = reverse @ Einv
print("Product in reverse order =\n", reverse)
print("\n ➜ The product in reverse order is still a unit lower triangular matrix, but it is not L (see the check in Step 4)")
# --- Step 4: check: product in order = L (assembled directly), product in reverse order ≠ L ---
print_step(4, "Check: product in order = the L assembled directly from the multipliers of each column; product in reverse order ≠ L")
L_direct = np.eye(n)
for k, Einv in enumerate(Einv_list):
L_direct[k+1:, k] = Einv[k+1:, k] # enter the elimination multipliers of each E_k⁻¹ directly
print("\nL assembled directly =\n", L_direct)
compare_print(
"Maximum difference ‖product in order − L_direct‖∞",
f"{np.max(np.abs(forward - L_direct)):.2e}",
"should be 0 (or machine precision): the product in order is L"
)
compare_print(
"Maximum difference ‖product in reverse order − L_direct‖∞",
f"{np.max(np.abs(reverse - L_direct)):.2e}",
"should be clearly nonzero: in the product in reverse order the multipliers contaminate one another, so it is not L"
)
assert np.allclose(forward, L_direct), "The product in order should equal the directly assembled L"
assert not np.allclose(reverse, L_direct), "The product in reverse order should not equal L"
print("\n ✓ Conclusion: L = E_0⁻¹E_1⁻¹⋯E_{n-1}⁻¹ (in order of increasing index); the product in reverse order is not.")
6.3.4 Complete Gaussian Elimination and the PLU Decomposition¶
As described in the previous subsection, the ideal LU decomposition is often not feasible in actual computation, because:
a diagonal entry may be zero and so cannot serve as the pivot entry
even when it is nonzero, a very small value leads to numerical instability
Practical Gaussian elimination therefore needs to introduce a row-exchange strategy (in engineering practice, column exchanges are also added to improve the stability of the algorithm further), which gives the PLU decomposition.
Permutation matrices and row exchanges
A permutation matrix is a matrix obtained by rearranging the rows (or the columns) of the identity matrix. Left-multiplying a matrix by a permutation matrix changes the order of the rows of that matrix.
An important property of permutation matrices is that their inverse equals their transpose: P−1=P⊤.
The early stages of the PLU decomposition are no different from the LU decomposition, until during the elimination we find that rows k and i (i>k) need to be interchanged. Because row exchanges and row eliminations are interleaved, the key question is: how can the PLU form of the decomposition be preserved while they are interleaved?
Block representation of the current state.
Suppose elimination has been completed for the first k rows; the state of the matrix can then be written as
The result is still a unit lower triangular matrix—only the positions of m10 and m20 have been swapped (corresponding to the exchange of row 1 and row 2).
The general case. The structure of Qk guarantees that it permutes only the row indices after step k and does not touch the completed k×k upper-left block L00; and the lower-right block is I before the row exchange and is still I afterward (rows and columns are interchanged together, which leaves the identity matrix unchanged). Therefore
Pk+1=QkPk: the updated accumulated permutation matrix
Lk′=QkLkQk: the updated unit lower triangular matrix
Uk′=QkUk: the partially upper triangular matrix after the row exchange, which is also partially in REF
Key insight. By this clever device we convert “row exchanges during the computation” into “row exchanges before the computation begins,” always maintaining the standard form PA=LU. In other words, if we knew the final permutation matrix P in advance, the ideal LU decomposition of PA would need no row exchanges at all!
Practical significance
The advantages of the PLU decomposition are:
Numerical stability: choosing the pivot entry of largest absolute value greatly reduces the accumulation of rounding errors
General applicability: it applies to every matrix and is not restricted by singularity
Computational efficiency: although row exchange operations are added, the overall complexity is still O(mnmin(m,n))
Theoretical completeness: it provides a reliable numerical foundation for the many applications of linear algebra
In practice, the PLU decomposition is the core tool for solving systems of linear equations, computing determinants, finding inverse matrices, and similar operations.
The following code implements the algorithm above; please read it alongside the steps
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ──────────────────────────────────────────────
def print_header(title):
print("=" * 60)
print(f" {title}")
print("=" * 60)
def print_step(step, desc):
print(f"\n▶ Step {step}: {desc}")
print("-" * 40)
def compare_print(label, actual, expected_desc):
print(f"\n[{label}]")
print(f" Computed value: {actual}")
print(f" Expected value: {expected_desc}")
# ──────────────────────────────────────────────
def flexible_plu_decomposition(A, tol=1e-12):
"""
PLU decomposition with flexible pivoting (pivots not restricted to diagonal positions)
Parameters
----------
A : input matrix, of shape m×n
tol : numerical tolerance for deciding whether a pivot is zero
Returns
-------
P : permutation matrix, of shape m×m
L : unit lower triangular matrix, of shape m×m
U : approximate REF matrix at this tolerance (a REF in exact arithmetic), of shape m×n
pivot_positions : list of pivot positions [(row, col), ...]
row_swaps : record of row exchanges [(row_i, row_j), ...]
success : Boolean indicating whether the decomposition succeeded
"""
m, n = A.shape
P = np.eye(m)
L = np.eye(m)
U = A.copy().astype(float)
row_swaps = []
pivot_positions = []
current_row = 0
current_col = 0
while current_row < m and current_col < n:
# --- Step 1: find the pivot of largest absolute value in the current column ---
col_segment = U[current_row:m, current_col]
rel_best = np.argmax(np.abs(col_segment))
best_row = current_row + rel_best
best_val = U[best_row, current_col]
# --- Step 2: if the pivot is too small, skip this column ---
if abs(best_val) < tol:
current_col += 1
continue
# --- Step 3: if the pivot is not in the current row, perform a row exchange ---
if best_row != current_row:
row_swaps.append((current_row, best_row))
U[[current_row, best_row]] = U[[best_row, current_row]]
P[[current_row, best_row]] = P[[best_row, current_row]]
if current_row > 0:
L[[current_row, best_row], :current_row] = \
L[[best_row, current_row], :current_row]
pivot_positions.append((current_row, current_col))
pivot = U[current_row, current_col]
# --- Step 4: eliminate, clearing the entries below the pivot ---
for i in range(current_row + 1, m):
if abs(U[i, current_col]) >= tol:
mult = U[i, current_col] / pivot
L[i, current_row] = mult
U[i, current_col:] = U[i, current_col:] - mult * U[current_row, current_col:]
U[i, current_col] = 0.0
current_row += 1
current_col += 1
return P, L, U, pivot_positions, row_swaps, True
# ──────────────────────────────────────────────
def verify_flexible_plu(A, P, L, U, tol=1e-10):
"""
Verify the PLU decomposition: PA = LU, equivalently A = P^T·LU
Returns bool: whether the verification passes
"""
err_PA_LU = np.max(np.abs(P @ A - L @ U))
err_A_rec = np.max(np.abs(A - P.T @ L @ U))
print(f" ‖PA − LU‖∞ = {err_PA_LU:.2e}")
print(f" ‖A − P^T LU‖∞ = {err_A_rec:.2e}")
return max(err_PA_LU, err_A_rec) < tol
def analyze_flexible_plu(P, L, U, pivot_positions, row_swaps):
"""
Analyze the result of the PLU decomposition: number of row exchanges, pivot positions, rank of the matrix
"""
m, n = U.shape
print(f" Number of row exchanges: {len(row_swaps)}")
if row_swaps:
for k, (r1, r2) in enumerate(row_swaps):
print(f" Exchange {k+1}: row {r1} ↔ row {r2}")
print(f" Numerical rank at this tolerance (number of pivots): {len(pivot_positions)}")
for k, (r, c) in enumerate(pivot_positions):
print(f" Pivot {k}: position ({r}, {c}), value = {U[r, c]:.4f}")
pivot_cols = [c for _, c in pivot_positions]
if len(pivot_cols) > 1:
ok = all(pivot_cols[i] < pivot_cols[i+1] for i in range(len(pivot_cols)-1))
print(f" Pivot column positions strictly increasing (REF property): {ok}")
# ──────────────────────────────────────────────
def run_case(label, A):
"""Run the complete procedure for a single test case"""
print(f"\n{'─'*60}")
print(f" {label}")
print(f"{'─'*60}")
print("A =\n", A)
P, L, U, positions, swaps, _ = flexible_plu_decomposition(A)
print("\nP =\n", P)
print("\nL =\n", L)
print("\nU =\n", U)
print("\n [Verification]")
passed = verify_flexible_plu(A, P, L, U)
print("\n [Analysis]")
analyze_flexible_plu(P, L, U, positions, swaps)
compare_print(
"Verification result",
"passed ✓" if passed else "failed ✗",
"PA = LU and A = P^T LU; the errors should be close to machine precision"
)
# ============================================================
print_header("Example 6.4 | PLU Decomposition with Flexible Pivoting")
# ============================================================
# --- Step 1: a square matrix that needs a row exchange ---
print_step(1, "A square matrix A₁ that needs a row exchange")
A1 = np.array([
[1.0, 1.0, 0.0, 0.0],
[1.0/2, 1.0/2, 1.0, 0.0],
[1.0/3, 4.0/3, 1.0, 0.0],
[1.0/4, 1.0/4, 1.0, 1.0],
])
run_case("Case 1: a square matrix that needs a row exchange", A1)
# --- Step 2: a matrix whose pivots are not on the diagonal ---
print_step(2, "A matrix A₂ whose pivots move off the diagonal after elimination")
A2 = np.array([
[1.0, 2.0, 0.0, 3.0],
[2.0, 4.0, 1.0, 7.0],
[1.0, 2.0, 2.0, 8.0],
])
run_case("Case 2: pivots off the diagonal (column 1 becomes a zero column after elimination)", A2)
# --- Step 3: a rectangular matrix (number of rows < number of columns) ---
print_step(3, "A rectangular matrix A₃ (3×4)")
A3 = np.array([
[1.0, 2.0, 3.0, 4.0],
[2.0, 1.0, 4.0, 3.0],
[1.0, 3.0, 2.0, 5.0],
])
run_case("Case 3: a rectangular matrix", A3)
# --- Step 4: a tiny pivot, which needs a row exchange to maintain numerical stability ---
print_step(4, "A matrix A₄ with a tiny pivot (demonstrating numerical stability)")
A4 = np.array([
[1e-10, 1.0, 2.0],
[1.0, 2.0, 1.0],
[2.0, 1.0, 3.0],
])
run_case("Case 4: tiny pivot → a row exchange ensures numerical stability", A4)
6.3.6 A Unified Framework for Solving Linear Systems and Inverting Matrices¶
The PLU decomposition obtained in the previous subsection can not only solve the system of linear equations Ax=b but also compute the inverse of a matrix efficiently. In fact, finding the inverse matrix is equivalent to solving n special systems of linear equations simultaneously, Axi=ei (where ei is the i-th column vector of the identity matrix), so that Gaussian elimination, the LU decomposition, and the PLU decomposition provide a unified theoretical and computational framework for solving systems of linear equations and computing inverse matrices.
Features of Gaussian elimination. In essence, it transforms the system of linear equations Ax=b by a sequence of row operations into an equivalent upper triangular system Ux=b′. The “advantage” of Gaussian elimination is that only the final simplified system needs to be recorded, not the transformation process, which makes it look “simple and intuitive”; even high school students can master it. However, because the transformation process is not recorded, every elimination step must be carried out again each time a new constant vector b comes along.
Features of the LU decomposition. This is an idealized method of matrix decomposition, which assumes that the pivots always lie on the diagonal during the elimination, so that no row exchanges are needed. Under this assumption we can express the elimination process as A=LU, where L records all the elimination multipliers and U is an upper triangular matrix. To solve a system of linear equations we only need to solve two triangular systems in turn: first Ly=b and then Ux=y. The advantage of this method is that the decomposition is computed once and can be used many times, while the constant vector b stays unchanged. However, the existence of the LU decomposition depends on the lucky condition that “no pivot is zero,” which does not always hold in practical applications.
Features of the PLU decomposition. This is the most complete and robust method of matrix decomposition: it introduces a permutation matrix P to record the necessary row exchanges and obtains PA=LU. For the system of linear equations Ax=b, we turn to solving PAx=Pb, and then solve Ly=Pb and Ux=y in turn. Since the permutation merely rearranges the entries of a vector, computing Pb is very efficient. The PLU decomposition combines the efficiency of decomposition methods with the numerical stability of row exchanges, and applies to almost every matrix.
Comparing the three methods. In terms of computational efficiency, for a single solve Gaussian elimination has a slight advantage because it avoids storing matrices; but when many systems with the same coefficient matrix must be solved, the advantage of the PLU decomposition is significant, because the decomposition is computed only once. In terms of numerical stability, Gaussian elimination and the PLU decomposition behave identically when they use the same pivoting strategy. In memory usage, Gaussian elimination has a slight advantage because it does not need to store the matrices L and P, but this difference is usually not decisive in modern computing environments.
Developments in modern numerical methods. Note that the three methods above all belong to the category of direct methods, which obtain the exact solution (within numerical precision) in finitely many steps. For large sparse systems, however, modern numerical linear algebra prefers iterative methods, such as the conjugate gradient method and GMRES. These methods need no explicit matrix decomposition; instead they approach the true solution by repeatedly improving an approximate solution, and they have better computational complexity and memory efficiency on large problems. The development of iterative methods, together with advances in preconditioning, has made them the mainstream approach to solving systems of linear equations in modern scientific computing.
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ──────────────────────────────────────────────
def print_header(title):
print("=" * 60)
print(f" {title}")
print("=" * 60)
def print_step(step, desc):
print(f"\n▶ Step {step}: {desc}")
print("-" * 40)
def compare_print(label, actual, expected_desc):
print(f"\n[{label}]")
print(f" Computed value: {actual}")
print(f" Expected value: {expected_desc}")
# ============================================================
print_header("Verification | Hand Solution of a System of Three Linear Equations (Compare with the Example Above)")
# ============================================================
# --- Step 1: define the coefficient matrix A and the right-hand-side vector b ---
print_step(1, "Define A and b")
A = np.array([
[ 2.0, 1.0, 1.0],
[-1.0, 1.0, -1.0],
[ 3.0, -1.0, -1.0]
])
b = np.array([8.0, -2.0, 5.0])
print("A =\n", A)
print("b =", b)
# --- Step 2: solve with np.linalg.solve ---
print_step(2, "Solve Ax = b with np.linalg.solve")
x = np.linalg.solve(A, b)
x_hand = np.array([13/5, 17/10, 11/10]) # hand-computed result x0=13/5, x1=17/10, x2=11/10
print("numerical solution x =", x)
print("hand solution x_hand =", x_hand)
# --- Step 3: comparison and residual ---
print_step(3, "Compare with the hand solution and check the residual")
compare_print("Maximum difference ‖x − x_hand‖∞", f"{np.max(np.abs(x - x_hand)):.2e}", "should be close to machine precision")
compare_print("Residual ‖Ax − b‖∞", f"{np.max(np.abs(A @ x_hand - b)):.2e}", "should be 0 when the hand solution is substituted back into the original system")
assert np.allclose(x, x_hand)
print("\n ✓ The hand-computed result x = (13/5, 17/10, 11/10) passes the check.")
Gaussian elimination can also be used to find the inverse of a matrix. For an n×n matrix A, construct the augmented matrix [A∣In], where In is the n×n identity matrix, and then apply Gaussian elimination to this augmented matrix. Note, however, that, as mentioned in the previous subsection, modern industry rarely computes inverse matrices; most of the time it prefers to solve the system of linear equations directly with iterative methods, and when an inverse matrix is really needed, the PLU decomposition is usually used, since it keeps the overall error under better control. Computing the inverse by Gaussian elimination in this subsection is therefore aimed more at understanding than at practical application.
6.4 Practical Applications of Systems of Linear Equations¶
Systems of linear equations have a wide range of applications in the real world; they appear everywhere from engineering design to economic analysis, and from network flows to chemical equilibrium.
6.4.1 Matrix Representation of Network Flow Problems¶
Network flow problems show the important value of systems of linear equations in practical applications. We will see that flow conservation in a network can be expressed as a system of linear equations Af=b, and that, because of the special properties of the incidence matrix, this system usually has infinitely many solutions. The essence of a network flow problem is optimization over the degrees of freedom in the space formed by these solutions (in particular, the kernel of the homogeneous system Af=0), searching for the optimal distribution of flow subject to the capacity constraints.
Matrix Representation of the Flow Conservation Equations¶
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
# Fonts: the English edition needs no CJK font, so nothing is downloaded when _lang == 'en';
# the Chinese editions use this same block to fetch Noto Sans TC/SC where no CJK font is installed
import os, urllib.request
import matplotlib.font_manager as fm
_lang = 'en'
_cjk = ['Microsoft JhengHei', 'PingFang TC', 'Noto Sans CJK TC', 'Noto Sans TC']
_have = {f.name for f in fm.fontManager.ttflist}
if _lang != 'en' and not _have & set(_cjk):
_font = os.path.join(os.path.expanduser('~'), '.cache', 'fonts', 'NotoSansTC.ttf')
try:
if not os.path.exists(_font):
os.makedirs(os.path.dirname(_font), exist_ok=True)
urllib.request.urlretrieve('https://github.com/google/fonts/raw/main/ofl/notosanstc/NotoSansTC%5Bwght%5D.ttf', _font + '.part')
os.replace(_font + '.part', _font)
fm.fontManager.addfont(_font)
_have.add('Noto Sans TC')
except OSError as err:
print('Could not download the CJK font; Chinese text in figures may not display:', err)
plt.rcParams['font.family'] = [f for f in _cjk if f in _have] + ['DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
# book color palette
C_TEAL, C_ORANGE, C_LIGHT_VIOLET = '#2AD2C9', '#FF5D47', '#AB82C5'
# Create the list of nodes, including the source 's' and the sink 't'
nodes = ['s', 'a', 'b', 'c', 't']
# Define the incidence matrix to represent the flow conservation equations
# rows: nodes ['s', 'a', 'b', 'c', 't']
# columns: edges [('s','a'), ('s','b'), ('a','b'), ('a','c'), ('b','c'), ('b','t'), ('c','t')]
# values: +1 means inflow (the node is the end of the edge), -1 means outflow (the node is the start of the edge), 0 means not connected
# —— exactly the same convention as def-incidence-matrix in the text
incidence_matrix = np.array([
[-1, -1, 0, 0, 0, 0, 0], # s: outflow on edges (s,a) and (s,b)
[ 1, 0, -1, -1, 0, 0, 0], # a: inflow (s,a), outflow (a,b) and (a,c)
[ 0, 1, 1, 0, -1, -1, 0], # b: inflow (s,b) and (a,b), outflow (b,c) and (b,t)
[ 0, 0, 0, 1, 1, 0, -1], # c: inflow (a,c) and (b,c), outflow (c,t)
[ 0, 0, 0, 0, 0, 1, 1] # t: inflow (b,t) and (c,t)
])
# Define the edges and their capacities
edges = [('s', 'a', 10), ('s', 'b', 8),
('a', 'b', 5), ('a', 'c', 4),
('b', 'c', 6), ('b', 't', 7),
('c', 't', 10)]
# Create the directed graph from the matrix
G = nx.DiGraph()
G.add_nodes_from(nodes)
for u, v, capacity in edges:
G.add_edge(u, v, capacity=capacity)
# Set the node positions for visualization
pos = {'s': (0, 1), 'a': (1, 2), 'b': (1, 0), 'c': (2, 1), 't': (3, 1)}
# Compute the maximum flow
max_flow_value, flow_dict = nx.maximum_flow(G, 's', 't')
# Set the edge labels in the format "flow/capacity"
edge_labels = {}
for u, v, capacity in edges:
flow = flow_dict[u][v]
edge_labels[(u, v)] = f"{flow}/{capacity}"
# Give the source and the sink special colors (book color palette)
node_colors = []
for node in nodes:
if node == 's':
node_colors.append(C_TEAL) # source: teal
elif node == 't':
node_colors.append(C_ORANGE) # sink: orange
else:
node_colors.append(C_LIGHT_VIOLET) # intermediate nodes: light violet
# Draw the network flow graph
plt.figure(figsize=(10, 6))
nx.draw(G, pos, with_labels=True, node_color=node_colors,
node_size=500, font_weight='bold', arrows=True)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.title("Network Flow Graph Built from the Flow Conservation Matrix")
plt.axis('off')
plt.tight_layout()
plt.show()
# Print the incidence matrix
print("Incidence matrix of flow conservation:")
print(incidence_matrix)
# Print the flow conservation equation of each node
print("\nFlow conservation equation of each node:")
edge_names = [f"f({u},{v})" for u, v, _ in edges]
for i, node in enumerate(nodes):
eq_parts = []
for j, (u, v, _) in enumerate(edges):
if incidence_matrix[i, j] == 1:
eq_parts.append(f"+{edge_names[j]}")
elif incidence_matrix[i, j] == -1:
eq_parts.append(f"-{edge_names[j]}")
equation = "".join(eq_parts)
if node == 's':
print(f"Source '{node}': {equation} = -v (v is the total flow; by convention the net flow at the source is negative)")
elif node == 't':
print(f"Sink '{node}': {equation} = v (v is the total flow; by convention the net flow at the sink is positive)")
else:
print(f"Node '{node}': {equation} = 0 (inflow = outflow)")
# The maximum flow was already computed in the plotting part
print(f"\nMaximum flow value: {max_flow_value}")
# Show the flow on each edge
print("\nOptimal flow assignment:")
flow_vector = []
for u, v, _ in edges:
flow = flow_dict[u][v]
flow_vector.append(flow)
print(f"f({u},{v}) = {flow}")
# Verify flow conservation: check A·f = b node by node
flow_vector = np.array(flow_vector)
node_net_flows = incidence_matrix @ flow_vector
# Construct the demand vector b by the def-incidence-matrix convention: source -v, sink +v, intermediate nodes 0
b = np.zeros(len(nodes))
b[nodes.index('s')] = -max_flow_value
b[nodes.index('t')] = max_flow_value
print("\nVerify flow conservation (A·f = b):")
for i, node in enumerate(nodes):
if node == 's':
print(f"Source '{node}': net flow = {node_net_flows[i]} (should be -{max_flow_value})")
elif node == 't':
print(f"Sink '{node}': net flow = {node_net_flows[i]} (should be {max_flow_value})")
else:
print(f"Node '{node}': net flow = {node_net_flows[i]} (should be 0)")
assert np.allclose(incidence_matrix @ flow_vector, b), "A·f should equal b"
print("\n ✓ A·f = b checked node by node: the code, the definition, and the example use the same sign convention (+1 = inflow).")
6.4.2 A Case Study in Balancing Chemical Equations¶
Balancing a chemical reaction equation is in essence solving a system of linear equations.
In mechanics problems, a body acted on by several forces remains in equilibrium, and this can be described by a system of linear equations.
Systems of linear equations have wide applications in every field. By turning network flow problems, chemical equilibrium, and mechanical equilibrium into systems of linear equations, we can analyze and solve these problems systematically. The theoretical framework and computational tools provided by linear algebra make complex problems solvable, showing the power of mathematics in practical applications.
This chapter started from the premise that “row operations are invertible linear transformations” and built a complete theoretical framework for solving systems of linear equations. §6.1 laid the foundation for the legitimacy of the operations: the three elementary row operations (row scaling, row exchange, row elimination) are each equivalent to left multiplication by an elementary matrix, and so they are invertible, change neither the solution set nor the rank of the matrix. The theorem on the preservation of subspace dimensions gives the most general statement—left or right multiplication by any invertible matrix preserves the dimensions of the row space and the column space—which guarantees the legitimacy of every subsequent reduction step.
§6.2 used this guarantee and, centered on the block structure [Ir0F0] of the RREF, gave a complete description of the geometric form of the solutions. The rank r determines all the dimensions of the four subspaces; the number of free variables, n−r, determines the dimension of the solution space; and the relation yp=−Fyf between the pivot variables and the free variables directly gives a parametric representation of the solution space. The three forms—no solution, a unique solution, infinitely many solutions—become clear at a glance from this block viewpoint, and the abstract question of “solvability” is reduced to a computable “comparison of ranks.”
§6.3 provided the algorithmic tools for bringing a general matrix to RREF. The LU decomposition is elegant in form and records the elimination multipliers directly as the lower triangular matrix L, but it depends on the precondition that no pivot is zero; the PLU decomposition records the row exchanges in a permutation matrix P, applies to every matrix, and the resulting U is itself already in row echelon form. The two share the same design philosophy: decompose once, use many times, so that when the same coefficient matrix meets different right-hand-side vectors, only two triangular systems need to be solved for each. §6.4 then used three very different applications—network flow, chemical balancing, and mechanical equilibrium—to show the breadth of Ax=b as a unified modeling language.
The link between this chapter and Chapter 5 (block matrices) lies in the block structure of the RREF itself: Chapter 5 established the algebraic rules for operating with block matrices, and this chapter applies them to reveal the rank structure of a matrix and the geometric form of the solution space. The rigorous derivation of the PLU decomposition relies on the multiplicative properties of block matrices, in particular on the key property that “the commutation relation between permutation matrices and lower triangular matrices allows the row exchanges to be moved, equivalently, to the very front of the decomposition”—a conclusion that becomes natural only from the block viewpoint.
The connection with Chapter 7 (determinants) is even more direct: the basic way of computing a determinant is precisely to bring the matrix to upper triangular form by row operations, and the effects of the row operations on the determinant—scaling multiplies it by c, an exchange changes its sign, elimination leaves its value unchanged—all follow naturally from the elementary-matrix framework of this chapter. The concept of rank in this chapter also finally merges with the theory of dimension of abstract linear spaces in Chapter 4: the concrete computational statement that “the number of pivots equals the dimension of the image” is exactly the dimension theorem of Chapter 4 expressed in the language of coordinates. Looking further ahead, the idea of the LU decomposition will reappear in Chapter 8 (the eigenvalue problem) and Chapter 11 (the singular value decomposition) in the forms of the Schur decomposition, the QR decomposition, and the SVD; “factoring a matrix into a product of factors of simpler structure” is a central theme that runs throughout the book.
Systems of linear equations are the oldest problem in linear algebra and also its most central one. The significance of this chapter lies not in supplying techniques for “solving equations,” but in revealing the answers to two essential questions: “what determines solvability?” and “what is the structure of the solutions?”; and in raising the arithmetic elimination already recorded in The Nine Chapters on the Mathematical Art to a rigorous algebraic language. Elimination operations are isomorphisms, the RREF is a canonical form, and the solution set is an affine subspace—these three sentences capture the full depth of this chapter. If you keep reading with this viewpoint, you will find that every later chapter answers the same deeper question: how to find a “canonical form” in which the essential structure of a matrix is laid bare. This chapter is the first complete demonstration of this theme.
