Matrix Inverses
Concept
A square matrix $\mathbf{A}$ is invertible if there exists another square matrix $\mathbf{A}^{-1}$, called its inverse, such that $\mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$. The inverse "undoes" the matrix: it is the linear algebra analogue of dividing by a number. Not every matrix has one — matrices without an inverse are called singular, and detecting them is as important as computing the inverse itself. The general algorithm for finding $\mathbf{A}^{-1}$ is Gauss-Jordan elimination: row-reduce the augmented matrix $[\mathbf{A} \mid \mathbf{I}]$ until the left half becomes $\mathbf{I}$, and the right half is then $\mathbf{A}^{-1}$.
Why This Matters
In ch05 you learned that a matrix is a transformation. The inverse is the transformation that undoes it — scale by 2, then scale by ½; rotate by θ, then rotate by −θ. This unlocks the most important use of matrices: solving linear systems. A system $\mathbf{A}\mathbf{x} = \mathbf{b}$ has the solution $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$, which is how least-squares regression (the normal equations $\mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b}$), covariance computations, and countless ML formulas get solved in practice. When a matrix has no inverse, the system it represents has no unique solution — knowing why is the first step toward understanding rank, determinants (ch07), and the geometry of underdetermined problems.
Mathematical Notation
The inverse
For a square matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$, the inverse is written $\mathbf{A}^{-1}$ (read "A inverse") and satisfies:
where $\mathbf{I}_n$ is the $n \times n$ identity matrix (introduced in ch03). The superscript $-1$ is not an exponent — $\mathbf{A}^{-1}$ never means $1/\mathbf{A}$, because there is no matrix division. The inverse is defined by this property alone.
Motivation: the defining equation $A A^{-1} = I$ is the target we verify in every test — if $\mathbf{A}\mathbf{A}^{-1} \approx \mathbf{I}$ within floating-point error, the inverse is correct.
Invertible vs. singular
A square matrix is invertible (also called non-singular) if its inverse exists. It is singular if no inverse exists. The key geometric fact: an invertible transformation maps distinct vectors to distinct vectors (it never collapses a direction), while a singular transformation squashes some direction to zero — information is destroyed, so it cannot be undone.
Motivation: every algorithm must handle the singular case explicitly. Our implementation returns None for singular matrices rather than panicking or producing garbage.
The augmented matrix
The augmented matrix $[\mathbf{A} \mid \mathbf{B}]$ is a single matrix that stacks $\mathbf{B}$ to the right of $\mathbf{A}$, separated by a vertical bar:
It is just a bookkeeping device: the two matrices share the same rows, so row operations applied to the whole augmented matrix affect both halves simultaneously.
Motivation: Gauss-Jordan computes $\mathbf{A}^{-1}$ by building $[\mathbf{A} \mid \mathbf{I}]$ and row-reducing until the left half is $\mathbf{I}$ — the right half is then forced to be $\mathbf{A}^{-1}$.
Elementary row operations
Three operations on rows never change the solution set of a system (equivalently, they never change invertibility):
- Row swap $R_a \leftrightarrow R_b$ — exchange two rows.
- Row scaling $R_a \leftarrow c \cdot R_a$ — multiply a row by a nonzero constant $c$.
- Row addition $R_a \leftarrow R_a + c \cdot R_b$ — add a multiple of one row to another.
Motivation: these three operations are exactly the private helper methods swap_rows, scale_row, and add_scaled_row in the implementation — the code is a direct translation of the math.
The 2×2 inverse formula
For a $2 \times 2$ matrix there is a closed form:
The denominator $ad - bc$ is the determinant — the subject of the next chapter. If it is zero, the matrix is singular and the formula fails, which is exactly when $\mathbf{A}^{-1}$ does not exist.
Intuition
Inverse as "undo"
A matrix transforms the whole space: every vector $\mathbf{v}$ becomes $\mathbf{A}\mathbf{v}$. If $\mathbf{A}$ is invertible, the transformation is a bijection — every output comes from exactly one input — so the process can be reversed. $\mathbf{A}^{-1}$ is that reversal: apply $\mathbf{A}$, then $\mathbf{A}^{-1}$, and you end up exactly where you started.
This matches the transformations you built in ch05:
- Scale by $(s_x, s_y)$ → inverse scales by $(1/s_x, 1/s_y)$.
- Rotate by $\theta$ → inverse rotates by $-\theta$. (And in fact for rotations, $\mathbf{R}^{-1} = \mathbf{R}^T$ — the transpose undoes the rotation. The test
test_inverse_rotation_is_transposechecks this.) - Reflect across the $x$-axis → inverse is the same reflection (reflecting twice returns you to the start).
Singular = collapse
A singular matrix collapses at least one direction: some nonzero vector $\mathbf{v}$ is mapped to $\mathbf{0}$. Once two different inputs map to the same output, no inverse can exist — you can't tell which input you started from. Think of it as a projection: a photo of a 3D scene onto a flat image can't be "un-projected" because the depth information was destroyed.
What Gauss-Jordan is really doing
The augmented matrix $[\mathbf{A} \mid \mathbf{I}]$ says: "I want to find $\mathbf{X}$ such that $\mathbf{A}\mathbf{X} = \mathbf{I}$." Applying a row operation to the whole augmented matrix is legal because row operations correspond to multiplying both sides by an invertible matrix on the left — they preserve the equation. When the left half reads $\mathbf{I}$, the right half must be the unique $\mathbf{X}$ that works, i.e. $\mathbf{A}^{-1}$.
Worked Examples
Example 1: 2×2 with the closed form
Verify with the defining property (this is test_inverse_2x2_known):
Example 2: 3×3 Gauss-Jordan trace
Take $\mathbf{A} = \begin{pmatrix} 1 & 2 & 3 \\ 0 & 1 & 4 \\ 5 & 6 & 0 \end{pmatrix}$. Write $[\mathbf{A} \mid \mathbf{I}]$ and row-reduce:
Column 0: the largest $|$entry$|$ in column 0 is $|5|$ (row 2), so swap row 0 ↔ row 2, then scale row 0 by $1/5$:
Column 1: pivot is already 1 in row 1. Eliminate above and below:
Column 2: scale row 2 by $-5$ to make the pivot 1, then eliminate above:
The left half is $\mathbf{I}$, so the right half is the inverse:
Notice this example swaps rows even though no pivot is exactly zero: at column 0, $|5| > |1|$, so partial pivoting brings the larger-magnitude entry to the pivot position. This is a numerical-stability choice — larger pivots keep floating-point error small (the true zero-pivot rescue appears in the 2×2 example $\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$). This exact trace is verified in test_inverse_general_3x3.
Example 3: a singular matrix
Row 2 is twice row 1, so Gauss-Jordan reaches a column where every candidate pivot is zero — the matrix collapses a direction and has no inverse. Our code returns None (test_singular_returns_none).
Rust Implementation
Add a new crate to your workspace:
cd code && cargo new --lib --edition 2024 ch06-matrix-inverses
This crate builds on the Matrix struct from ch03. Open code/ch06-matrix-inverses/src/lib.rs and start by copying the full Matrix impl from ch03-linear-algebra-matrices/src/lib.rs (as you did in ch04). Then add the following methods and free functions:
/// A matrix is stored in row-major order as a flat Vec<f64>.
///
/// Row-major means the first row occupies indices 0..cols,
/// the second row occupies indices cols..2*cols, and so on.
/// The element at (row, col) is at index row * cols + col.
#[derive(Debug, Clone, PartialEq)]
pub struct Matrix {
data: Vec<f64>,
rows: usize,
cols: usize,
}
impl Matrix {
/// Create a new matrix from a flat Vec<f64> in row-major order.
///
/// # Panics
/// Panics if data.len() != rows * cols.
pub fn new(data: Vec<f64>, rows: usize, cols: usize) -> Self {
assert_eq!(
data.len(),
rows * cols,
"Data length {} does not match dimensions {}×{}",
data.len(),
rows,
cols
);
Matrix { data, rows, cols }
}
/// Number of rows.
pub fn rows(&self) -> usize {
self.rows
}
/// Number of columns.
pub fn cols(&self) -> usize {
self.cols
}
/// Shape as (rows, cols).
pub fn shape(&self) -> (usize, usize) {
(self.rows, self.cols)
}
/// Access the element at (row, col) — 0-indexed.
///
/// # Panics
/// Panics if row or col is out of bounds.
pub fn get(&self, row: usize, col: usize) -> f64 {
assert!(row < self.rows, "Row {} out of bounds (rows={})", row, self.rows);
assert!(col < self.cols, "Col {} out of bounds (cols={})", col, self.cols);
self.data[row * self.cols + col]
}
/// Matrix multiplication: self * other.
///
/// For A ∈ ℝ^{m×n} and B ∈ ℝ^{n×p}, the result C ∈ ℝ^{m×p}.
/// C_{ij} = Σ_{k=0}^{n-1} A_{ik} * B_{kj}
///
/// # Panics
/// Panics if self.cols != other.rows (incompatible dimensions).
pub fn multiply(&self, other: &Matrix) -> Matrix {
assert_eq!(
self.cols, other.rows,
"Incompatible dimensions: {}×{} vs {}×{}",
self.rows, self.cols, other.rows, other.cols
);
let m = self.rows;
let n = self.cols;
let p = other.cols;
let mut data = Vec::with_capacity(m * p);
for i in 0..m {
for j in 0..p {
let mut sum = 0.0;
for k in 0..n {
sum += self.get(i, k) * other.get(k, j);
}
data.push(sum);
}
}
Matrix { data, rows: m, cols: p }
}
/// Transpose the matrix: flip rows and columns.
///
/// If A ∈ ℝ^{m×n}, then A^T ∈ ℝ^{n×m} with (A^T)_{ij} = A_{ji}.
pub fn transpose(&self) -> Matrix {
let mut data = Vec::with_capacity(self.rows * self.cols);
for j in 0..self.cols {
for i in 0..self.rows {
data.push(self.get(i, j));
}
}
Matrix {
data,
rows: self.cols,
cols: self.rows,
}
}
/// Create an n×n identity matrix.
pub fn identity(n: usize) -> Matrix {
let mut data = vec![0.0; n * n];
for i in 0..n {
data[i * n + i] = 1.0;
}
Matrix {
data,
rows: n,
cols: n,
}
}
/// Apply this matrix to a vector (as a column matrix).
///
/// If this matrix is m×n, the input vector must have n components,
/// and the result has m components.
///
/// Mathematically: T(v) = A * v where A is this matrix.
pub fn transform(&self, v: &[f64]) -> Vec<f64> {
assert_eq!(
self.cols,
v.len(),
"Matrix cols {} doesn't match vector length {}",
self.cols,
v.len()
);
let mut result = vec![0.0; self.rows];
for i in 0..self.rows {
let mut sum = 0.0;
for j in 0..self.cols {
sum += self.get(i, j) * v[j];
}
result[i] = sum;
}
result
}
/// ── Matrix inverses ──
/// Build the augmented matrix [self | rhs] by appending the columns of
/// `rhs` to the right of `self`.
///
/// For A ∈ ℝ^{m×n} and B ∈ ℝ^{m×p}, the result is the m×(n+p) matrix
/// whose first n columns are A and whose last p columns are B.
/// Writing the two side by side: [A | B].
///
/// This is the working form used by Gauss-Jordan elimination, where we
/// reduce [A | I] to [I | A⁻¹].
///
/// # Panics
/// Panics if self.rows != rhs.rows (must share the same row count).
pub fn augmented(&self, rhs: &Matrix) -> Matrix {
assert_eq!(
self.rows, rhs.rows,
"Augmented matrices must have the same row count ({} vs {})",
self.rows, rhs.rows
);
let mut data = Vec::with_capacity(self.rows * (self.cols + rhs.cols));
for i in 0..self.rows {
for j in 0..self.cols {
data.push(self.get(i, j));
}
for j in 0..rhs.cols {
data.push(rhs.get(i, j));
}
}
Matrix {
data,
rows: self.rows,
cols: self.cols + rhs.cols,
}
}
/// Compute the inverse of this square matrix via Gauss-Jordan
/// elimination with partial pivoting.
///
/// Returns `None` if the matrix is singular (a zero pivot appears),
/// since a singular matrix has no inverse.
///
/// # Panics
/// Panics if the matrix is not square.
pub fn inverse(&self) -> Option<Matrix> {
assert_eq!(
self.rows, self.cols,
"Inverse requires a square matrix (got {}×{})",
self.rows, self.cols
);
let n = self.rows;
let mut aug = self.augmented(&Matrix::identity(n));
for col in 0..n {
// Partial pivoting: among rows col..n, find the one with the
// largest |value| in this column and swap it into the pivot row.
// This keeps pivots away from zero and improves numerical
// stability (avoids dividing by tiny numbers).
let mut pivot = col;
for row in (col + 1)..n {
if aug.get(row, col).abs() > aug.get(pivot, col).abs() {
pivot = row;
}
}
// If every candidate pivot is ~0, the column is all zeros below
// the diagonal, so the matrix is singular. The 1e-10 threshold is
// an absolute one: it assumes a well-scaled matrix (entries of
// order 1), which is the convention used throughout this course.
// A matrix with all entries below ~1e-10 would be falsely
// reported singular — acceptable for a teaching implementation.
if aug.get(pivot, col).abs() < 1e-10 {
return None;
}
aug.swap_rows(col, pivot);
// Normalize the pivot row: divide through by the pivot so the
// pivot entry becomes 1.
let pivot_value = aug.get(col, col);
aug.scale_row(col, 1.0 / pivot_value);
// Eliminate this column from every other row, leaving a 0 above
// and below the pivot.
for row in 0..n {
if row != col {
let factor = aug.get(row, col);
aug.add_scaled_row(row, col, -factor);
}
}
}
// The left half is now I and the right half is A⁻¹.
let mut data = Vec::with_capacity(n * n);
for i in 0..n {
for j in 0..n {
data.push(aug.get(i, n + j));
}
}
Some(Matrix {
data,
rows: n,
cols: n,
})
}
/// ── Elementary row operations (private helpers) ──
/// Swap rows `a` and `b` in place.
fn swap_rows(&mut self, a: usize, b: usize) {
if a == b {
return;
}
for j in 0..self.cols {
let tmp = self.data[a * self.cols + j];
self.data[a * self.cols + j] = self.data[b * self.cols + j];
self.data[b * self.cols + j] = tmp;
}
}
/// Multiply row `row` by `factor` in place.
fn scale_row(&mut self, row: usize, factor: f64) {
for j in 0..self.cols {
self.data[row * self.cols + j] *= factor;
}
}
/// Add `factor` times row `src` to row `dst` in place:
/// row_dst ← row_dst + factor * row_src.
fn add_scaled_row(&mut self, dst: usize, src: usize, factor: f64) {
for j in 0..self.cols {
self.data[dst * self.cols + j] += factor * self.data[src * self.cols + j];
}
}
}
#[cfg(test)]
mod tests {
use super::*;
fn approx_eq(a: f64, b: f64) -> bool {
(a - b).abs() < 1e-10
}
fn assert_matrix_eq(a: &Matrix, b: &Matrix) {
assert_eq!(a.shape(), b.shape());
for i in 0..a.rows() {
for j in 0..a.cols() {
assert!(
approx_eq(a.get(i, j), b.get(i, j)),
"Mismatch at ({}, {}): {} vs {}",
i,
j,
a.get(i, j),
b.get(i, j)
);
}
}
}
fn assert_is_identity(m: &Matrix) {
for i in 0..m.rows() {
for j in 0..m.cols() {
let expected = if i == j { 1.0 } else { 0.0 };
assert!(
approx_eq(m.get(i, j), expected),
"Not identity at ({}, {}): {}",
i,
j,
m.get(i, j)
);
}
}
}
#[test]
fn test_augmented_shape_and_values() {
// [A | I] for a 2×2: columns 0..2 are A, columns 2..4 are I.
let a = Matrix::new(vec![4.0, 7.0, 2.0, 6.0], 2, 2);
let aug = a.augmented(&Matrix::identity(2));
assert_eq!(aug.shape(), (2, 4));
assert!(approx_eq(aug.get(0, 0), 4.0));
assert!(approx_eq(aug.get(0, 1), 7.0));
assert!(approx_eq(aug.get(1, 0), 2.0));
assert!(approx_eq(aug.get(1, 1), 6.0));
// right half is the identity
assert!(approx_eq(aug.get(0, 2), 1.0));
assert!(approx_eq(aug.get(0, 3), 0.0));
assert!(approx_eq(aug.get(1, 2), 0.0));
assert!(approx_eq(aug.get(1, 3), 1.0));
}
#[test]
#[should_panic(expected = "same row count")]
fn test_augmented_row_mismatch_panics() {
let a = Matrix::new(vec![1.0, 2.0, 3.0], 1, 3);
let b = Matrix::new(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
a.augmented(&b);
}
#[test]
fn test_inverse_2x2_known() {
// The classic example: A = [[4, 7], [2, 6]].
// det = 4*6 - 7*2 = 10, so
// A⁻¹ = 1/10 * [[6, -7], [-2, 4]] = [[0.6, -0.7], [-0.2, 0.4]].
let a = Matrix::new(vec![4.0, 7.0, 2.0, 6.0], 2, 2);
let inv = a.inverse().expect("A should be invertible");
let expected = Matrix::new(vec![0.6, -0.7, -0.2, 0.4], 2, 2);
assert_matrix_eq(&inv, &expected);
}
#[test]
fn test_inverse_2x2_known_with_zero_pivot() {
// A has a 0 in the (0,0) position — without pivoting, Gauss-Jordan
// would divide by zero. Partial pivoting swaps rows to fix this.
// A = [[0, 1], [1, 0]] is its own inverse.
let a = Matrix::new(vec![0.0, 1.0, 1.0, 0.0], 2, 2);
let inv = a.inverse().expect("A should be invertible");
let expected = Matrix::new(vec![0.0, 1.0, 1.0, 0.0], 2, 2);
assert_matrix_eq(&inv, &expected);
}
#[test]
fn test_inverse_diagonal() {
// The inverse of a diagonal matrix is diagonal with reciprocal entries.
let a = Matrix::new(vec![2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 4.0], 3, 3);
let inv = a.inverse().expect("diagonal with nonzero entries is invertible");
let expected = Matrix::new(
vec![0.5, 0.0, 0.0, 0.0, 1.0 / 3.0, 0.0, 0.0, 0.0, 0.25],
3,
3,
);
assert_matrix_eq(&inv, &expected);
}
#[test]
fn test_inverse_rotation_is_transpose() {
// From ch05: for a rotation R(θ), R⁻¹ = Rᵀ (rotating back by -θ).
// Verify that Gauss-Jordan produces exactly the transpose.
let theta: f64 = 0.7;
let r = Matrix::new(
vec![
theta.cos(),
-theta.sin(),
theta.sin(),
theta.cos(),
],
2,
2,
);
let inv = r.inverse().expect("rotation is invertible");
assert_matrix_eq(&inv, &r.transpose());
}
#[test]
fn test_inverse_general_3x3() {
// A = [[1, 2, 3], [0, 1, 4], [5, 6, 0]] — a generic invertible 3×3.
// Known result (hand-computed Gauss-Jordan trace, requires a row swap):
// A⁻¹ = [[-24, 18, 5], [20, -15, -4], [-5, 4, 1]]
let a = Matrix::new(vec![1.0, 2.0, 3.0, 0.0, 1.0, 4.0, 5.0, 6.0, 0.0], 3, 3);
let inv = a.inverse().expect("A should be invertible");
let expected = Matrix::new(
vec![-24.0, 18.0, 5.0, 20.0, -15.0, -4.0, -5.0, 4.0, 1.0],
3,
3,
);
assert_matrix_eq(&inv, &expected);
// The defining property A * A⁻¹ = I and A⁻¹ * A = I.
assert_is_identity(&a.multiply(&inv));
assert_is_identity(&inv.multiply(&a));
}
#[test]
fn test_inverse_of_inverse() {
// (A⁻¹)⁻¹ = A.
let a = Matrix::new(vec![1.0, 2.0, 3.0, 0.0, 1.0, 4.0, 5.0, 6.0, 0.0], 3, 3);
let inv = a.inverse().expect("A should be invertible");
let back = inv.inverse().expect("A⁻¹ should also be invertible");
assert_matrix_eq(&back, &a);
}
#[test]
fn test_inverse_identity_is_identity() {
let i = Matrix::identity(4);
let inv = i.inverse().expect("identity is invertible");
assert_matrix_eq(&inv, &i);
}
#[test]
fn test_solve_linear_system() {
// Solve Ax = b using x = A⁻¹b.
// A = [[2, 1], [1, 3]], b = [5, 10].
// det = 2*3 - 1*1 = 5, so A⁻¹ = 1/5 [[3, -1], [-1, 2]].
// x = A⁻¹b = 1/5 [3*5 - 10, -5 + 2*10] = 1/5 [5, 15] = [1, 3].
let a = Matrix::new(vec![2.0, 1.0, 1.0, 3.0], 2, 2);
let b = vec![5.0, 10.0];
let inv = a.inverse().expect("A should be invertible");
let x = inv.transform(&b);
assert!(approx_eq(x[0], 1.0), "x[0] = {}", x[0]);
assert!(approx_eq(x[1], 3.0), "x[1] = {}", x[1]);
// Sanity: A * x should recover b.
let ax = a.transform(&x);
assert!(approx_eq(ax[0], 5.0));
assert!(approx_eq(ax[1], 10.0));
}
#[test]
fn test_singular_returns_none() {
// Rows are linearly dependent: row2 = 2 * row1.
let a = Matrix::new(vec![1.0, 2.0, 2.0, 4.0], 2, 2);
assert!(a.inverse().is_none());
}
#[test]
fn test_zero_matrix_returns_none() {
let a = Matrix::new(vec![0.0, 0.0, 0.0, 0.0], 2, 2);
assert!(a.inverse().is_none());
}
#[test]
#[should_panic(expected = "square matrix")]
fn test_inverse_non_square_panics() {
let a = Matrix::new(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 2, 3);
let _ = a.inverse();
}
}
Run the tests:
cargo test -p ch06-matrix-inverses
You should see:
running 13 tests
test tests::test_augmented_row_mismatch_panics - should panic ... ok
test tests::test_augmented_shape_and_values ... ok
test tests::test_inverse_2x2_known ... ok
test tests::test_inverse_2x2_known_with_zero_pivot ... ok
test tests::test_inverse_diagonal ... ok
test tests::test_inverse_general_3x3 ... ok
test tests::test_inverse_identity_is_identity ... ok
test tests::test_inverse_non_square_panics - should panic ... ok
test tests::test_inverse_of_inverse ... ok
test tests::test_inverse_rotation_is_transpose ... ok
test tests::test_singular_returns_none ... ok
test tests::test_solve_linear_system ... ok
test tests::test_zero_matrix_returns_none ... ok
test result: ok. 13 passed; 0 failed; 0 ignored; 0 measured; 0 filtered
Walkthrough
-
augmented(&self, rhs)— Builds $[\mathbf{A} \mid \mathbf{B}]$ by copying each row of $\mathbf{A}$ followed by the same row of $\mathbf{B}$. The row counts must match (both matrices share rows). This is pure bookkeeping — no arithmetic, just layout. -
inverse(&self) -> Option<Matrix>— The heart of the chapter. It builds $[\mathbf{A} \mid \mathbf{I}]$ and runs one pass per column:- Find the pivot — the row at or below the diagonal with the largest $|$value$|$ in this column.
- Check for singularity — if every candidate pivot is ~0, the column is zeroed out below the diagonal, meaning the matrix collapses a direction. Return
None. - Swap, normalize, eliminate — swap the pivot row into position, divide the row by the pivot (making the pivot 1), then subtract multiples of the pivot row from every other row to zero out the column.
When every column is done, the left half is $\mathbf{I}$ and the right half is $\mathbf{A}^{-1}$. Why
Option? Because singularity is a legitimate outcome, not an error —Noneis the honest answer for "this matrix has no inverse." -
Partial pivoting — The step most textbook descriptions skip, and the reason the algorithm actually works in practice. Without it, a perfectly invertible matrix like $\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$ would crash on a zero pivot (
test_inverse_2x2_known_with_zero_pivot). Picking the largest-magnitude pivot also keeps floating-point error small: dividing by a tiny number amplifies rounding noise. -
The three private helpers —
swap_rows,scale_row,add_scaled_roware exactly the three elementary row operations from the notation section, written as in-place&mut selfmethods. Keeping them private hides the mechanics;inversereads as a clean description of the algorithm. -
transform(&self, v)— Carried over from ch05. It lets us use the inverse: solve $\mathbf{A}\mathbf{x} = \mathbf{b}$ withx = a.inverse().unwrap().transform(&b), as intest_solve_linear_system. -
The 2×2 formula is a preview —
1/(ad−bc)is the determinant. For larger matrices there's no simple closed form, which is exactly why we need the general Gauss-Jordan algorithm. ch07 will formalize the determinant and give a direct singularity test.
Verification
| Test | What it checks | Mathematical invariant |
|---|---|---|
test_augmented_shape_and_values | $[\mathbf{A} \mid \mathbf{I}]$ layout | Augmented matrix construction |
test_augmented_row_mismatch_panics | Row-count check on augmented | Both halves share rows |
test_inverse_2x2_known | Known 2×2 inverse via the closed form | $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$ |
test_inverse_2x2_known_with_zero_pivot | Zero pivot rescued by row swap | Partial pivoting correctness |
test_inverse_diagonal | Diagonal matrix inverse | diag reciprocal entries |
test_inverse_general_3x3 | Hand-computed 3×3 inverse + both sides | $\mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ |
test_inverse_of_inverse | Inverting twice returns the original | $(\mathbf{A}^{-1})^{-1} = \mathbf{A}$ |
test_inverse_rotation_is_transpose | Rotation inverse equals its transpose | $\mathbf{R}(\theta)^{-1} = \mathbf{R}(\theta)^T$ |
test_inverse_identity_is_identity | Identity inverts to itself | $\mathbf{I}^{-1} = \mathbf{I}$ |
test_solve_linear_system | Solves $\mathbf{A}\mathbf{x} = \mathbf{b}$ via $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$ | $\mathbf{A}\mathbf{A}^{-1}\mathbf{b} = \mathbf{b}$ |
test_singular_returns_none | Linearly dependent rows → no inverse | Singular detection |
test_zero_matrix_returns_none | Zero matrix → no inverse | Singular detection |
test_inverse_non_square_panics | Non-square input panics | Inverses require square matrices |
Key Takeaways
- A square matrix $\mathbf{A}$ is invertible when a matrix $\mathbf{A}^{-1}$ exists with $\mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$; the inverse undoes the transformation (ch05).
- Gauss-Jordan elimination computes $\mathbf{A}^{-1}$ by row-reducing $[\mathbf{A} \mid \mathbf{I}]$ to $[\mathbf{I} \mid \mathbf{A}^{-1}]$ using three elementary row operations: swap, scale, add.
- Partial pivoting (swapping in the largest-magnitude pivot) is what makes the algorithm robust — it avoids zero pivots and reduces floating-point error.
- A matrix with no inverse is singular; the algorithm detects this when every candidate pivot in a column is ~0 and returns
None. - The inverse solves linear systems: $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$ solves $\mathbf{A}\mathbf{x} = \mathbf{b}$ — the foundation for least squares and much of numerical ML.
- The 2×2 closed form reveals the determinant ($ad - bc$) as the singularity test — the subject of the next chapter.