MIT 18.06SC Lecture 2: Elimination with Matrices

Linear Algebra

MIT 18.06SC

Gaussian Elimination

Author

Chao Ma

Published

October 1, 2025

This recap of MIT 18.06SC Lecture 2 explores Gaussian elimination—the systematic algorithm that transforms any linear system into an easily solvable upper triangular form.

📓 For a deeper dive with additional exercises and analysis, see the complete notebook on GitHub.

From Linear Systems to Upper Triangular Form

How do we actually solve a system of linear equations like this?

\[\begin{cases} 2u + v + w = 5 \\ 4u - 6v = -2 \\ -2u + 7v + 2w = 9 \end{cases}\]

The answer is Gaussian elimination—a systematic process that transforms the coefficient matrix into an upper triangular form, where solutions can be read off directly through back substitution. This fundamental algorithm underpins much of numerical linear algebra and is essential for understanding how computers solve linear systems.

Quick Reference: Understanding Elimination

For context on the mathematical foundations of Gaussian elimination, see the Lecture 2 summary.

The Two-Step Process:

Forward Elimination: Transform $A$ into upper triangular $U$ using row operations
- For each pivot position, eliminate all entries below it
- Row operation: $\text{Row}_i \leftarrow \text{Row}_i - m_{ij} \cdot \text{Row}_j$ where $m_{ij} = a_{ij}/\text{pivot}$
Back Substitution: Solve $Ux = c$ from bottom to top
- Start with the last equation (only one unknown)
- Work upward, substituting known values

Key insight: Each elimination step can be represented as multiplication by an elimination matrix $E_{ij}$, connecting row operations to matrix multiplication.

🔬 Exercise: Implement Gaussian Elimination from Scratch

Let’s implement the complete Gaussian elimination algorithm and verify it works correctly.

Show code

import numpy as np
from IPython.display import display, Markdown

# Set random seed for reproducibility
np.random.seed(42)

def matrix_to_latex(matrix, precision=2):
    """Convert numpy matrix to LaTeX format."""
    if len(matrix.shape) == 1:
        # Vector
        elements = " \\\\ ".join([f"{x:.{precision}f}" for x in matrix])
        return f"$$\\begin{{bmatrix}}{elements}\\end{{bmatrix}}$$"
    else:
        # Matrix
        rows = []
        for row in matrix:
            row_str = " & ".join([f"{x:.{precision}f}" for x in row])
            rows.append(row_str)
        matrix_str = " \\\\ ".join(rows)
        return f"$$\\begin{{bmatrix}}{matrix_str}\\end{{bmatrix}}$$"

print("✓ Setup complete")

✓ Setup complete

Algorithm Overview

The algorithm consists of two phases:

Forward Elimination: Transform $A$ into upper triangular form $U$
- For each column $j$ from 0 to $n-1$:
  - Pivot = $A[j, j]$
  - For each row $i$ below pivot ($i > j$):
    - Compute multiplier: $m_{ij} = A[i, j] / \text{pivot}$
    - Row operation: $A[i, :] = A[i, :] - m_{ij} \cdot A[j, :]$
    - Update RHS: $b[i] = b[i] - m_{ij} \cdot b[j]$
Back Substitution: Solve $Ux = c$ from bottom to top
- Start from last row: $x[n-1] = c[n-1] / U[n-1, n-1]$
- For each row $i$ from $n-2$ down to $0$:
  - $x[i] = (c[i] - \sum_{j=i+1}^{n-1} U[i,j] \cdot x[j]) / U[i,i]$

Implementation

Show code

def gaussian_elimination(A, b):
    """
    Solve the linear system Ax = b using Gaussian elimination.

    Parameters:
    -----------
    A : numpy.ndarray, shape (n, n)
        Coefficient matrix
    b : numpy.ndarray, shape (n,)
        Right-hand side vector

    Returns:
    --------
    x : numpy.ndarray, shape (n,)
        Solution vector
    """
    A_copy = A.copy()
    b_copy = b.copy()

    rows, cols = A.shape
    assert rows == cols, "Matrix must be square"

    for i in range(rows - 1):
        multiplier = A_copy[i + 1:, i] / A_copy[i, i]
        for index, m in enumerate(multiplier):
            A_copy[i + 1 + index] = A_copy[i + 1 + index] - m * A_copy[i]
            assert abs(A_copy[i + 1 + index, i]) < 1e-10, "Upper triangular matrix expected"
            b_copy[i + 1 + index] = b_copy[i + 1 + index] - m * b_copy[i]

    U = A_copy
    print("Upper triangular matrix U:")
    display(Markdown(matrix_to_latex(U)))

    # Now our matrix is upper triangular
    # Solve for x from bottom to top
    x = np.zeros(rows)
    for i in range(rows - 1, -1, -1):
        x[i] = (b_copy[i] - np.dot(A_copy[i][i + 1:], x[i + 1:])) / A_copy[i][i]

    return x

print("✓ Function defined")

✓ Function defined

Test on Lecture Example

Test on the system from the lecture:

\[\begin{cases} 2u + v + w = 5 \\ 4u - 6v = -2 \\ -2u + 7v + 2w = 9 \end{cases}\]

Expected solution: $u = 1, v = 1, w = 2$

Show code

# Define the system from lecture
A = np.array([
    [2, 1, 1],
    [4, -6, 0],
    [-2, 7, 2]
], dtype=float)

b = np.array([5, -2, 9], dtype=float)

print("System to solve:")
print("A =")
print(A)
print("\nb =")
print(b)

# Solve using our implementation
x_our = gaussian_elimination(A, b)

# Solve using NumPy
x_numpy = np.linalg.solve(A, b)

# Compare results
print("\n" + "=" * 60)
print("Results:")
print("=" * 60)
print(f"Our implementation:      x = {x_our}")
print(f"NumPy (np.linalg.solve): x = {x_numpy}")
print("=" * 60)

# Verify solution
residual_our = np.linalg.norm(A @ x_our - b)
residual_numpy = np.linalg.norm(A @ x_numpy - b)

print(f"\nResidual (our):   ||Ax - b|| = {residual_our:.2e}")
print(f"Residual (NumPy): ||Ax - b|| = {residual_numpy:.2e}")

# Check difference
difference = np.linalg.norm(x_our - x_numpy)
print(f"\nDifference: ||x_our - x_numpy|| = {difference:.2e}")

if difference < 1e-10:
    print("\n✅ Solutions match! Implementation is correct.")
else:
    print("\n❌ Solutions differ. Check implementation.")

System to solve:
A =
[[ 2.  1.  1.]
 [ 4. -6.  0.]
 [-2.  7.  2.]]

b =
[ 5. -2.  9.]
Upper triangular matrix U:

\[\begin{bmatrix}2.00 & 1.00 & 1.00 \\ 0.00 & -8.00 & -2.00 \\ 0.00 & 0.00 & 1.00\end{bmatrix}\]


============================================================
Results:
============================================================
Our implementation:      x = [1. 1. 2.]
NumPy (np.linalg.solve): x = [1. 1. 2.]
============================================================

Residual (our):   ||Ax - b|| = 0.00e+00
Residual (NumPy): ||Ax - b|| = 0.00e+00

Difference: ||x_our - x_numpy|| = 0.00e+00

✅ Solutions match! Implementation is correct.

Visualize Upper Triangular Matrix

Let’s test on a larger 10×10 system to see the upper triangular structure more clearly.

Show code

# Create a 10x10 random system
np.random.seed(123)
A_large = np.random.randn(10, 10)
b_large = np.random.randn(10)

print("Original 10×10 matrix A:")
display(Markdown(matrix_to_latex(A_large)))

print("\n" + "=" * 70)

# Solve using our implementation
x_large = gaussian_elimination(A_large, b_large)

print("=" * 70)
print("\n✓ Upper triangular structure achieved!")

Original 10×10 matrix A:

\[\begin{bmatrix}-1.09 & 1.00 & 0.28 & -1.51 & -0.58 & 1.65 & -2.43 & -0.43 & 1.27 & -0.87 \\ -0.68 & -0.09 & 1.49 & -0.64 & -0.44 & -0.43 & 2.21 & 2.19 & 1.00 & 0.39 \\ 0.74 & 1.49 & -0.94 & 1.18 & -1.25 & -0.64 & 0.91 & -1.43 & -0.14 & -0.86 \\ -0.26 & -2.80 & -1.77 & -0.70 & 0.93 & -0.17 & 0.00 & 0.69 & -0.88 & 0.28 \\ -0.81 & -1.73 & -0.39 & 0.57 & 0.34 & -0.01 & 2.39 & 0.41 & 0.98 & 2.24 \\ -1.29 & -1.04 & 1.74 & -0.80 & 0.03 & 1.07 & 0.89 & 1.75 & 1.50 & 1.07 \\ -0.77 & 0.79 & 0.31 & -1.33 & 1.42 & 0.81 & 0.05 & -0.23 & -1.20 & 0.20 \\ 0.47 & -0.83 & 1.16 & -1.10 & -2.12 & 1.04 & -0.40 & -0.13 & -0.84 & -1.61 \\ 1.26 & -0.69 & 1.66 & 0.81 & -0.31 & -1.09 & -0.73 & -1.21 & 2.09 & 0.16 \\ 1.15 & -1.27 & 0.18 & 1.18 & -0.34 & 1.03 & -1.08 & -1.36 & 0.38 & -0.38\end{bmatrix}\]


======================================================================
Upper triangular matrix U:

\[\begin{bmatrix}-1.09 & 1.00 & 0.28 & -1.51 & -0.58 & 1.65 & -2.43 & -0.43 & 1.27 & -0.87 \\ -0.00 & -0.72 & 1.31 & 0.30 & -0.08 & -1.47 & 3.72 & 2.46 & 0.21 & 0.93 \\ -0.00 & 0.00 & 3.22 & 1.07 & -1.89 & -3.94 & 10.50 & 5.69 & 1.36 & 1.35 \\ -0.00 & 0.00 & 0.00 & 0.82 & -2.93 & -3.41 & 8.91 & 3.46 & 1.04 & -0.34 \\ 0.00 & 0.00 & 0.00 & 0.00 & 6.41 & 7.26 & -17.36 & -8.55 & -1.51 & 2.79 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.15 & -1.79 & -0.91 & -0.09 & -0.48 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 3.20 & 0.39 & -1.78 & -0.72 \\ -0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & -2.27 & 0.17 & 2.74 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.73 & -2.90 \\ -0.00 & -0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 3.24\end{bmatrix}\]

======================================================================

✓ Upper triangular structure achieved!

This implementation demonstrates the fundamental algorithm for solving linear systems, connecting the geometric view from Lecture 1 with the algebraic machinery of row operations and matrix elimination.

--- title: "MIT 18.06SC Lecture 2: Elimination with Matrices" author: "Chao Ma" date: "2025-10-01" categories: ["Linear Algebra", "MIT 18.06SC", "Gaussian Elimination"] code-fold: true code-summary: "Show code" --- *This recap of MIT 18.06SC Lecture 2 explores Gaussian elimination—the systematic algorithm that transforms any linear system into an easily solvable upper triangular form.* {{< video https://www.youtube.com/watch?v=QVKj3LADCnA >}} 📓 **For a deeper dive with additional exercises and analysis**, see the [complete notebook on GitHub](https://github.com/ickma2311/foundations/blob/main/MIT18.06SC/lecture02/exercises.ipynb). ## From Linear Systems to Upper Triangular Form How do we actually solve a system of linear equations like this? $$\begin{cases} 2u + v + w = 5 \\ 4u - 6v = -2 \\ -2u + 7v + 2w = 9 \end{cases}$$ The answer is **Gaussian elimination**—a systematic process that transforms the coefficient matrix into an upper triangular form, where solutions can be read off directly through back substitution. This fundamental algorithm underpins much of numerical linear algebra and is essential for understanding how computers solve linear systems. ### Quick Reference: Understanding Elimination For context on the mathematical foundations of Gaussian elimination, see the [Lecture 2 summary](https://github.com/ickma2311/foundations/blob/main/MIT18.06SC/lecture02/lecture02_elimination.md). **The Two-Step Process**: 1. **Forward Elimination**: Transform $A$ into upper triangular $U$ using row operations - For each pivot position, eliminate all entries below it - Row operation: $\text{Row}_i \leftarrow \text{Row}_i - m_{ij} \cdot \text{Row}_j$ where $m_{ij} = a_{ij}/\text{pivot}$ 2. **Back Substitution**: Solve $Ux = c$ from bottom to top - Start with the last equation (only one unknown) - Work upward, substituting known values **Key insight**: Each elimination step can be represented as multiplication by an **elimination matrix** $E_{ij}$, connecting row operations to matrix multiplication. ## 🔬 Exercise: Implement Gaussian Elimination from Scratch Let's implement the complete Gaussian elimination algorithm and verify it works correctly. ```{python} import numpy as np from IPython.display import display, Markdown # Set random seed for reproducibility np.random.seed(42) def matrix_to_latex(matrix, precision=2): """Convert numpy matrix to LaTeX format.""" if len(matrix.shape) == 1: # Vector elements = " \\\\ ".join([f"{x:.{precision}f}" for x in matrix]) return f"$$\\begin{{bmatrix}}{elements}\\end{{bmatrix}}$$" else: # Matrix rows = [] for row in matrix: row_str = " & ".join([f"{x:.{precision}f}" for x in row]) rows.append(row_str) matrix_str = " \\\\ ".join(rows) return f"$$\\begin{{bmatrix}}{matrix_str}\\end{{bmatrix}}$$" print("✓ Setup complete") ``` ### Algorithm Overview The algorithm consists of two phases: 1. **Forward Elimination**: Transform $A$ into upper triangular form $U$ - For each column $j$ from 0 to $n-1$: - Pivot = $A[j, j]$ - For each row $i$ below pivot ($i > j$): - Compute multiplier: $m_{ij} = A[i, j] / \text{pivot}$ - Row operation: $A[i, :] = A[i, :] - m_{ij} \cdot A[j, :]$ - Update RHS: $b[i] = b[i] - m_{ij} \cdot b[j]$ 2. **Back Substitution**: Solve $Ux = c$ from bottom to top - Start from last row: $x[n-1] = c[n-1] / U[n-1, n-1]$ - For each row $i$ from $n-2$ down to $0$: - $x[i] = (c[i] - \sum_{j=i+1}^{n-1} U[i,j] \cdot x[j]) / U[i,i]$ ### Implementation ```{python} def gaussian_elimination(A, b): """ Solve the linear system Ax = b using Gaussian elimination. Parameters: ----------- A : numpy.ndarray, shape (n, n) Coefficient matrix b : numpy.ndarray, shape (n,) Right-hand side vector Returns: -------- x : numpy.ndarray, shape (n,) Solution vector """ A_copy = A.copy() b_copy = b.copy() rows, cols = A.shape assert rows == cols, "Matrix must be square" for i in range(rows - 1): multiplier = A_copy[i + 1:, i] / A_copy[i, i] for index, m in enumerate(multiplier): A_copy[i + 1 + index] = A_copy[i + 1 + index] - m * A_copy[i] assert abs(A_copy[i + 1 + index, i]) < 1e-10, "Upper triangular matrix expected" b_copy[i + 1 + index] = b_copy[i + 1 + index] - m * b_copy[i] U = A_copy print("Upper triangular matrix U:") display(Markdown(matrix_to_latex(U))) # Now our matrix is upper triangular # Solve for x from bottom to top x = np.zeros(rows) for i in range(rows - 1, -1, -1): x[i] = (b_copy[i] - np.dot(A_copy[i][i + 1:], x[i + 1:])) / A_copy[i][i] return x print("✓ Function defined") ``` ### Test on Lecture Example Test on the system from the lecture: $$\begin{cases} 2u + v + w = 5 \\ 4u - 6v = -2 \\ -2u + 7v + 2w = 9 \end{cases}$$ Expected solution: $u = 1, v = 1, w = 2$ ```{python} # Define the system from lecture A = np.array([ [2, 1, 1], [4, -6, 0], [-2, 7, 2] ], dtype=float) b = np.array([5, -2, 9], dtype=float) print("System to solve:") print("A =") print(A) print("\nb =") print(b) # Solve using our implementation x_our = gaussian_elimination(A, b) # Solve using NumPy x_numpy = np.linalg.solve(A, b) # Compare results print("\n" + "=" * 60) print("Results:") print("=" * 60) print(f"Our implementation: x = {x_our}") print(f"NumPy (np.linalg.solve): x = {x_numpy}") print("=" * 60) # Verify solution residual_our = np.linalg.norm(A @ x_our - b) residual_numpy = np.linalg.norm(A @ x_numpy - b) print(f"\nResidual (our): ||Ax - b|| = {residual_our:.2e}") print(f"Residual (NumPy): ||Ax - b|| = {residual_numpy:.2e}") # Check difference difference = np.linalg.norm(x_our - x_numpy) print(f"\nDifference: ||x_our - x_numpy|| = {difference:.2e}") if difference < 1e-10: print("\n✅ Solutions match! Implementation is correct.") else: print("\n❌ Solutions differ. Check implementation.") ``` ### Visualize Upper Triangular Matrix Let's test on a larger 10×10 system to see the upper triangular structure more clearly. ```{python} # Create a 10x10 random system np.random.seed(123) A_large = np.random.randn(10, 10) b_large = np.random.randn(10) print("Original 10×10 matrix A:") display(Markdown(matrix_to_latex(A_large))) print("\n" + "=" * 70) # Solve using our implementation x_large = gaussian_elimination(A_large, b_large) print("=" * 70) print("\n✓ Upper triangular structure achieved!") ``` --- *This implementation demonstrates the fundamental algorithm for solving linear systems, connecting the geometric view from Lecture 1 with the algebraic machinery of row operations and matrix elimination.*