4. Matrix Algebra

Chapter 3 used matrices as compact bookkeeping for linear systems. This chapter treats them as algebraic objects in their own right: matrices can be added, scaled, multiplied, transposed, and sometimes inverted, and these operations obey rules that mostly, but not entirely, match the arithmetic of numbers. Those rules are the working language of the rest of the subject, from determinants and linear transformations to the factorizations that numerical libraries use to solve systems with millions of unknowns.

<a id="41-matrices-and-basic-operations"></a>

4.1 Matrices and Basic Operations

This subchapter defines matrices and their three basic operations (addition, scalar multiplication, and matrix multiplication), and identifies which laws of ordinary arithmetic survive and which fail.

Matrix Definition and Notation

Why it matters. Rectangular tables of numbers are everywhere: the coefficients of a linear system, the pixel intensities of a grayscale image, the distances between cities, the weights of a neural-network layer, the transition probabilities of a Markov chain. Treating a table as a single object with its own name and operations lets you reason about all of its entries at once, the same way vector notation lets you reason about all coordinates of a point at once.

Definition. An $m \times n$ matrix over a field $\mathbb{F}$ (usually $\mathbb{R}$ or $\mathbb{C}$) is a rectangular array of $mn$ elements of $\mathbb{F}$, called entries, arranged in $m$ rows and $n$ columns:

$$A = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \ a_{21} & a_{22} & \cdots & a_{2n} \ \vdots & \vdots & \ddots & \vdots \ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix} = [a_{ij}].$$

The pair $m \times n$ is the size (or shape) of $A$, read "$m$ by $n$". As in Chapter 3, the convention is always rows first, then columns.

Notation	Meaning
$a_{ij}$, $(A){ij}$, or $A{ij}$	Entry in row $i$, column $j$
$\mathbb{R}^{m \times n}$ or $M_{m \times n}(\mathbb{F})$	Set of all $m \times n$ matrices with real (or $\mathbb{F}$) entries
$M_n(\mathbb{F})$	Set of all square $n \times n$ matrices
$\operatorname{row}_i(A)$	Row $i$, a $1 \times n$ matrix
$\operatorname{col}_j(A)$	Column $j$, an $m \times 1$ matrix, i.e. a vector in $\mathbb{R}^m$
$A = [,\mathbf{c}_1\ \ \mathbf{c}_2\ \ \cdots\ \ \mathbf{c}_n,]$	$A$ written by its columns

Shapes with special names.

• A square matrix has $m = n$. Its main diagonal consists of the entries $a_{11}, a_{22}, \dots, a_{nn}$.
• A column vector is an $n \times 1$ matrix. From now on, a vector $\mathbf{x} \in \mathbb{R}^n$ is identified with the column vector of its components.
• A row vector is a $1 \times n$ matrix.
• A $1 \times 1$ matrix $[a]$ is identified with the scalar $a$.

Equality. Two matrices are equal if they have the same size and the same entry in every position. $\begin{bmatrix} 1 & 2 \end{bmatrix}$ and $\begin{bmatrix} 1 \ 2 \end{bmatrix}$ hold the same numbers but are different matrices, because their sizes differ.

Example. A tiny $3 \times 4$ grayscale image is a matrix whose entry $a_{ij}$ is the brightness of the pixel in row $i$, column $j$ (0 is black, 255 is white):

                 col 1  col 2  col 3  col 4
              [    0     50    100    150 ]   <- row 1
   A    =     [   50    100    150    200 ]   <- row 2
              [  100    150    200    255 ]   <- row 3
                                  ^
                                  a_23 = 150: row 2, column 3

   size 3 x 4;   row_2(A) = [50 100 150 200];   col_4(A) = (150, 200, 255)

Matrices can also be defined by a formula for their entries. The $3 \times 3$ matrix with $a_{ij} = i + j$ is $\begin{bmatrix} 2 & 3 & 4 \ 3 & 4 & 5 \ 4 & 5 & 6 \end{bmatrix}$.

import numpy as np

A = np.array([[  0,  50, 100, 150],
              [ 50, 100, 150, 200],
              [100, 150, 200, 255]])

print(A.shape)          # (3, 4): rows first, then columns
print(A[1, 2])          # 150: the entry a_23 (NumPy indices start at 0)
print(A[1, :])          # [ 50 100 150 200]: row 2
print(A[:, 3])          # [150 200 255]: column 4, returned as a flat 1-D array

# A matrix defined by a formula for its entries: a_ij = i + j, with 1-based i and j
B = np.fromfunction(lambda i, j: (i + 1) + (j + 1), (3, 3), dtype=int)
print(B)                # [[2 3 4] [3 4 5] [4 5 6]]

# Equality requires equal shapes as well as equal entries
row = np.array([[1, 2]])                   # shape (1, 2)
col = np.array([[1], [2]])                 # shape (2, 1)
print(row.shape == col.shape)              # False: same numbers, different matrices
print(np.array_equal(row, col))            # False

Key takeaways

• An $m \times n$ matrix has $m$ rows and $n$ columns; $a_{ij}$ is the entry in row $i$, column $j$.
• Vectors in $\mathbb{R}^n$ are treated as $n \times 1$ column matrices.
• Two matrices are equal only if they have the same size and agree entry by entry.
• Mathematical indices start at 1; NumPy indices start at 0, so $a_{ij}$ is A[i-1, j-1].

🧪 Practice

1. For $A = \begin{bmatrix} 4 & -1 & 0 & 7 \ 2 & 5 & 3 & -6 \end{bmatrix}$, give the size, $a_{14}$, $a_{23}$, $\operatorname{row}_2(A)$, and $\operatorname{col}_3(A)$.
2. Write out the $3 \times 4$ matrix with entries $a_{ij} = i - j$. Which entries are zero, and where do they lie?
3. Find all values of $x$ and $y$ for which $\begin{bmatrix} x + y & 2 \ 3 & x - y \end{bmatrix} = \begin{bmatrix} 5 & 2 \ 3 & 1 \end{bmatrix}$.
4. Interview: A $10{,}000 \times 10{,}000$ matrix is stored in memory as one long array. NumPy stores it row by row by default, and Fortran (and therefore LAPACK) stores it column by column. Why can summing each column be several times slower in one layout than the other? Hint: think about which entries are adjacent in memory and how CPU caches load memory in contiguous chunks.

Matrix Addition and Scalar Multiplication

Why it matters. Many natural operations on tables are performed entry by entry: adding this month's sales table to last month's, averaging two images, converting a whole table of prices from dollars to euros. These two operations make matrices behave exactly like vectors with $mn$ components, so everything learned about vector arithmetic in Section 2.2 carries over without change.

Definitions. Let $A = [a_{ij}]$ and $B = [b_{ij}]$ be matrices of the same size $m \times n$, and let $c$ be a scalar.

• The sum $A + B$ is the $m \times n$ matrix with entries $(A + B){ij} = a{ij} + b_{ij}$.
• The scalar multiple $cA$ is the $m \times n$ matrix with entries $(cA){ij} = c,a{ij}$.
• The negative is $-A = (-1)A$, and the difference is $A - B = A + (-B)$.

The sum of matrices of different sizes is undefined. There is no sensible way to add a $2 \times 3$ table to a $3 \times 2$ table.

Example.

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

• 2\begin{bmatrix} 2 & 1 & -1 \ 0 & -2 & 3 \end{bmatrix} = \begin{bmatrix} 1 & -2 & 0 \ 3 & 4 & 5 \end{bmatrix} + \begin{bmatrix} 4 & 2 & -2 \ 0 & -4 & 6 \end{bmatrix} = \begin{bmatrix} 5 & 0 & -2 \ 3 & 0 & 11 \end{bmatrix}. $$

Algebraic properties. Because each operation acts on every entry separately, each property below follows from the same property of numbers, applied entry by entry. For $m \times n$ matrices $A, B, C$ and scalars $c, d$:

Property	Statement
Commutativity	$A + B = B + A$
Associativity	$(A + B) + C = A + (B + C)$
Zero matrix	$A + O = A$, where $O$ is the $m \times n$ matrix of zeros
Additive inverse	$A + (-A) = O$
Distributivity over matrix sums	$c(A + B) = cA + cB$
Distributivity over scalar sums	$(c + d)A = cA + dA$
Compatibility of scalars	$c(dA) = (cd)A$
Unit scalar	$1A = A$

These are exactly the vector space axioms (Section 6.1), so $\mathbb{R}^{m \times n}$ is a vector space of dimension $mn$. Think of an $m \times n$ matrix as a vector with $mn$ entries that happens to be displayed in a rectangle.

Linear combinations of matrices. Combining the two operations gives linear combinations such as $(1 - t)A + tB$. For two images $A$ and $B$ of the same size, this is a cross-fade: $t = 0$ shows $A$, $t = 1$ shows $B$, and intermediate values blend them.

import numpy as np

A = np.array([[1, -2, 0], [3, 4, 5]])
B = np.array([[2, 1, -1], [0, -2, 3]])
print(A + 2 * B)                       # [[ 5  0 -2] [ 3  0 11]]

# Cross-fade between two 2 x 2 "images"
dark  = np.array([[ 0.,  0.], [  0.,   0.]])
light = np.array([[200., 100.], [100., 200.]])
for t in [0.0, 0.25, 1.0]:
    print(t, ((1 - t) * dark + t * light).tolist())   # 0.25 -> [[50.0, 25.0], [25.0, 50.0]]

# Mathematically, sizes must match. NumPy refuses incompatible shapes ...
try:
    A + np.ones((3, 2))
except ValueError as e:
    print("error:", e)                 # operands could not be broadcast together ...

# ... but silently "broadcasts" some mismatches, which is NOT matrix addition
v = np.array([10, 20, 30])             # a 1-D array of length 3
print(A + v)                           # [[11 18 30] [13 24 35]]: v was added to every row

The last line is a common source of bugs. NumPy's broadcasting stretches the shorter array to fit, so an expression that is mathematically undefined can run without error. When the shapes in your formulas matter, check them with assert A.shape == B.shape.

Key takeaways

• Addition and scalar multiplication act entry by entry; addition requires matrices of the same size.
• They satisfy the same laws as vector addition and scaling, so $\mathbb{R}^{m \times n}$ is a vector space of dimension $mn$.
• The zero matrix $O$ is the additive identity, and $-A$ is the additive inverse.
• NumPy broadcasting can silently produce results for mathematically undefined sums; check shapes explicitly.

🧪 Practice

1. Compute $3A - 2B$ for $A = \begin{bmatrix} 1 & 0 \ -1 & 2 \end{bmatrix}$ and $B = \begin{bmatrix} 2 & -3 \ 4 & 1 \end{bmatrix}$.
2. Solve for the matrix $X$: $2X + \begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix} = \begin{bmatrix} 5 & 0 \ -1 & 6 \end{bmatrix}$. Which properties from the table did you use?
3. Show that every $2 \times 2$ matrix can be written uniquely as a linear combination of $E_{11} = \begin{bmatrix} 1 & 0 \ 0 & 0 \end{bmatrix}$, $E_{12} = \begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$, $E_{21} = \begin{bmatrix} 0 & 0 \ 1 & 0 \end{bmatrix}$, $E_{22} = \begin{bmatrix} 0 & 0 \ 0 & 1 \end{bmatrix}$. What does this say about the dimension of $\mathbb{R}^{2 \times 2}$?
4. Interview: In a code review, you see loss = (predictions - targets) ** 2 where predictions has shape (64, 1) and targets has shape (64,). The code runs and the loss is suspiciously large. What happened? Hint: work out the shape of the difference under NumPy's broadcasting rules.

Matrix Multiplication

Why it matters. Matrix multiplication is not defined entry by entry, and the reason is that it has a more useful job: it represents doing one linear operation after another. Suppose a bakery's orders $\mathbf{x}$ (loaves of bread, cakes) determine the ingredients it needs, $\mathbf{y} = B\mathbf{x}$ (flour, sugar, eggs), and the ingredients determine the total cost and calories, $\mathbf{z} = A\mathbf{y}$. Then $\mathbf{z} = A(B\mathbf{x})$. Matrix multiplication is defined so that a single matrix $AB$ goes straight from orders to cost and calories: $(AB)\mathbf{x} = A(B\mathbf{x})$.

Definition. Let $A$ be $m \times n$ and $B$ be $n \times p$. The product $AB$ is the $m \times p$ matrix with entries

$$(AB){ij} = a{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{in}b_{nj} = \sum_{k=1}^{n} a_{ik}b_{kj}.$$

In words, entry $(i, j)$ of $AB$ is the dot product of row $i$ of $A$ with column $j$ of $B$. The product is defined only when the inner sizes match: the number of columns of $A$ must equal the number of rows of $B$.

           A            B            AB
        (m x n)      (n x p)      (m x p)
              \________/
           inner sizes must match;  the outer sizes give the shape of the result

   entry (i, j) of AB:

        row i of A  -->  [ a_i1  a_i2  ...  a_in ]       b_1j   <-- column j of B
                                                         b_2j
                         a_i1 b_1j + a_i2 b_2j + ... + a_in b_nj
                                                         b_nj

The matrix-vector product $A\mathbf{x}$ from Section 3.2 is the special case $p = 1$.

Example. With $A = \begin{bmatrix} 1 & 2 & 0 \ 0 & 1 & -1 \end{bmatrix}$ ($2 \times 3$) and $B = \begin{bmatrix} 3 & 1 \ 1 & 0 \ 2 & 4 \end{bmatrix}$ ($3 \times 2$):

$$AB = \begin{bmatrix} 1\cdot3 + 2\cdot1 + 0\cdot2 & 1\cdot1 + 2\cdot0 + 0\cdot4 \ 0\cdot3 + 1\cdot1 + (-1)\cdot2 & 0\cdot1 + 1\cdot0 + (-1)\cdot4 \end{bmatrix} = \begin{bmatrix} 5 & 1 \ -1 & -4 \end{bmatrix}.$$

The product $BA$ is also defined, but it is $3 \times 3$, so it cannot equal $AB$.

Four ways to see the same product. Each view is useful in different situations, and fluency in switching between them is one of the most valuable skills in linear algebra.

View	Formula	When it helps
Entry (row times column)	$(AB)_{ij} = \operatorname{row}_i(A) \cdot \operatorname{col}_j(B)$	Computing by hand, one entry at a time
Column	$AB = [,A\mathbf{b}_1\ \ A\mathbf{b}_2\ \ \cdots\ \ A\mathbf{b}_p,]$	Each column of $AB$ is $A$ applied to a column of $B$; columns of $AB$ are combinations of columns of $A$
Row	$\operatorname{row}_i(AB) = \operatorname{row}_i(A),B$	Rows of $AB$ are combinations of rows of $B$; used for row operations (Section 4.4)
Sum of outer products	$AB = \sum_{k=1}^{n} \operatorname{col}_k(A),\operatorname{row}_k(B)$	Low-rank structure, the SVD (Chapter 11), blocked algorithms

In the outer-product view, each term $\operatorname{col}_k(A),\operatorname{row}_k(B)$ is an $(m \times 1)(1 \times p) = m \times p$ matrix, and $AB$ is the sum of $n$ such matrices.

Cost. Each of the $mp$ entries needs $n$ multiplications, so the product costs $mnp$ multiplications. Two $1000 \times 1000$ matrices take $10^9$ multiplications. Optimized libraries (BLAS) reorder the loops to use the CPU cache efficiently and run many times faster than a naive triple loop, though the operation count is the same.

import numpy as np

A = np.array([[1, 2, 0], [0, 1, -1]])          # 2 x 3
B = np.array([[3, 1], [1, 0], [2, 4]])         # 3 x 2

def matmul(A, B):
    """Textbook definition: (AB)_ij = sum over k of a_ik * b_kj."""
    m, n = A.shape
    n2, p = B.shape
    assert n == n2, "inner sizes must match"
    C = np.zeros((m, p), dtype=A.dtype)
    for i in range(m):
        for j in range(p):
            for k in range(n):
                C[i, j] += A[i, k] * B[k, j]
    return C

print(matmul(A, B))                             # [[ 5  1] [-1 -4]]
print(A @ B)                                    # same: @ is matrix multiplication
print((B @ A).shape)                            # (3, 3): BA is defined but has a different size

# Column view: column j of AB is A times column j of B
print(np.array_equal(np.column_stack([A @ B[:, j] for j in range(2)]), A @ B))   # True
# Row view: row i of AB is row i of A times B
print(np.array_equal(np.vstack([A[i, :] @ B for i in range(2)]), A @ B))         # True
# Outer-product view: sum of (column k of A)(row k of B)
print(np.array_equal(sum(np.outer(A[:, k], B[k, :]) for k in range(3)), A @ B)) # True

# Pitfall: * is entrywise multiplication in NumPy, not the matrix product
S = np.array([[1, 2], [3, 4]])
print(S * S)                                    # [[ 1  4] [ 9 16]]
print(S @ S)                                    # [[ 7 10] [15 22]]

# Composition: orders -> ingredients -> (cost, calories)
Bk = np.array([[0.50, 0.30],                    # kg flour per loaf, per cake
               [0.02, 0.20],                    # kg sugar
               [0.00, 4.00]])                   # eggs
Ak = np.array([[1.20, 2.00, 0.25],              # dollars per kg flour, kg sugar, egg
               [3640, 3870, 70]])               # kcal per kg flour, kg sugar, egg
x = np.array([10, 3])                           # 10 loaves, 3 cakes
print(np.allclose(Ak @ (Bk @ x), (Ak @ Bk) @ x))   # True: one matrix does both steps
print((Ak @ Bk).round(2).tolist())              # [[0.64, 1.76], [1897.4, 2146.0]]: $ and kcal per loaf, per cake

Key takeaways

• $AB$ is defined when $A$ is $m \times n$ and $B$ is $n \times p$; the result is $m \times p$.
• Entry $(i, j)$ of $AB$ is row $i$ of $A$ dotted with column $j$ of $B$.
• The column, row, and outer-product views describe the same product and are often more illuminating than the entry formula.
• Matrix multiplication represents composition: $(AB)\mathbf{x} = A(B\mathbf{x})$.
• In NumPy, @ is the matrix product and * is the entrywise product.

🧪 Practice

1. Compute $AB$ and $BA$ for $A = \begin{bmatrix} 1 & -1 \ 2 & 0 \end{bmatrix}$ and $B = \begin{bmatrix} 3 & 1 \ 1 & 2 \end{bmatrix}$.
2. Let $A$ be $3 \times 5$, $B$ be $5 \times 2$, and $C$ be $2 \times 3$. Which of $AB$, $BA$, $BC$, $CA$, $ABC$, $CAB$, $A^2$ are defined, and what are their sizes?
3. Using the column view, explain why, if column 3 of $B$ is the sum of columns 1 and 2 of $B$, then column 3 of $AB$ is the sum of columns 1 and 2 of $AB$.
4. Interview: You need to compute $ABx$ where $A$ and $B$ are $n \times n$ and $x$ is a vector, with $n = 10{,}000$. Does it matter whether you compute $(AB)x$ or $A(Bx)$? Estimate both costs. Hint: count multiplications for a matrix-matrix product and for a matrix-vector product, and remember that both orders give the same answer.

Properties and Non-Commutativity

Why it matters. Manipulating matrix expressions is like manipulating algebraic expressions with numbers, except that a few familiar rules are false. Knowing exactly which rules survive keeps derivations correct; most errors in matrix algebra come from silently using a rule that only holds for numbers.

Laws that hold. For matrices of sizes that make each side defined, and any scalar $c$:

Law	Statement
Associativity	$(AB)C = A(BC)$
Left distributivity	$A(B + C) = AB + AC$
Right distributivity	$(A + B)C = AC + BC$
Scalars move freely	$c(AB) = (cA)B = A(cB)$
Identity	$I_m A = A = A I_n$ for $A$ of size $m \times n$ (Section 4.2)

Associativity means a product such as $ABCD$ needs no parentheses. It is not obvious from the entry formula, but the composition view makes it clear: applying $C$, then $B$, then $A$ is the same process however you group the steps. In indices, both sides of $(AB)C = A(BC)$ have entry $(i, j)$ equal to $\sum_k \sum_l a_{ik}b_{kl}c_{lj}$, with the two sums taken in different orders.

Laws that fail.

Rule for numbers	Fails for matrices	Counterexample or reason
$ab = ba$	$AB \ne BA$ in general	Shown below
$ab = 0 \Rightarrow a = 0$ or $b = 0$	Nonzero matrices can have a zero product	$\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}\begin{bmatrix} 1 & 0 \ 0 & 0 \end{bmatrix} = O$
$ab = ac,\ a \ne 0 \Rightarrow b = c$	Cancellation fails	If $AB = O$ with $B \ne O$, then $AB = A,O$
$(a + b)^2 = a^2 + 2ab + b^2$	$(A + B)^2 = A^2 + AB + BA + B^2$	Equal only when $AB = BA$ (Section 1.3)

Non-commutativity. $AB$ and $BA$ can differ in three increasingly subtle ways:

1. $BA$ may be undefined. If $A$ is $2 \times 3$ and $B$ is $3 \times 4$, then $AB$ is $2 \times 4$ but $BA$ does not exist.
2. Both may be defined but have different sizes, as in the $2 \times 3$ times $3 \times 2$ example of the previous topic.
3. Both may be square of the same size and still differ.

Example. Let $R = \begin{bmatrix} 0 & -1 \ 1 & 0 \end{bmatrix}$ rotate the plane by $90^\circ$ counterclockwise and $F = \begin{bmatrix} 1 & 0 \ 0 & -1 \end{bmatrix}$ reflect it across the $x$-axis. Then

$$RF = \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}, \qquad FR = \begin{bmatrix} 0 & -1 \ -1 & 0 \end{bmatrix}.$$

$RF$ means "reflect first, then rotate", because $RF\mathbf{x} = R(F\mathbf{x})$. It is the reflection across the line $y = x$. $FR$ means "rotate first, then reflect", and it is the reflection across $y = -x$. Following $\mathbf{e}_1 = (1, 0)$:

   RF e1:  reflect (1,0) -> (1,0),   then rotate -> (0, 1)
   FR e1:  rotate  (1,0) -> (0,1),   then reflect -> (0,-1)

              y
              |
       (0,1)  *  <- RF e1
              |
      --------+--------*------ x
              |      (1,0) = e1
              |
       (0,-1) *  <- FR e1

Getting dressed is the everyday analogy: socks then shoes is not the same as shoes then socks.

When matrices do commute. Some pairs commute: any square $A$ commutes with itself, with its powers, with $I$, and with every scalar matrix $cI$. Diagonal matrices of the same size commute with each other. When $AB = BA$, the binomial theorem and other identities from ordinary algebra are valid for $A$ and $B$. The commutator $[A, B] = AB - BA$ measures the failure to commute. It is zero exactly when $A$ and $B$ commute, and it is central in quantum mechanics, where non-commuting observables cannot be measured simultaneously with arbitrary precision.

import numpy as np

R = np.array([[0, -1], [1, 0]])      # rotate 90 degrees counterclockwise
F = np.array([[1, 0], [0, -1]])      # reflect across the x-axis
e1 = np.array([1, 0])

print(R @ F)                         # [[0 1] [1 0]]: reflection across y = x
print(F @ R)                         # [[ 0 -1] [-1  0]]: reflection across y = -x
print((R @ F) @ e1, (F @ R) @ e1)    # [0 1] [ 0 -1]

# Zero divisors: a product of nonzero matrices can be zero, in one order but not the other
A = np.array([[0, 1], [0, 0]])
B = np.array([[1, 0], [0, 0]])
print(A @ B)                         # [[0 0] [0 0]]
print(B @ A)                         # [[0 1] [0 0]]

# Cancellation fails: AC = AD although C != D
C = np.array([[1, 2], [3, 4]])
D = np.array([[7, 8], [3, 4]])       # differs from C only in row 1, which A ignores
print(np.array_equal(A @ C, A @ D), np.array_equal(C, D))    # True False

# Associativity and distributivity always hold (up to rounding for floats)
rng = np.random.default_rng(0)
X, Y, Z = (rng.standard_normal((3, 3)) for _ in range(3))
print(np.allclose((X @ Y) @ Z, X @ (Y @ Z)))      # True
print(np.allclose(X @ (Y + Z), X @ Y + X @ Z))    # True
print(np.allclose(X @ Y, Y @ X))                  # False: random matrices almost never commute

Key takeaways

• Matrix multiplication is associative and distributive, and scalars can be moved anywhere in a product.
• It is not commutative: $AB$ and $BA$ may be undefined, of different sizes, or simply different.
• $AB = O$ does not imply $A = O$ or $B = O$, and $AB = AC$ does not imply $B = C$.
• Identities from ordinary algebra such as $(A + B)^2 = A^2 + 2AB + B^2$ hold only when the matrices commute.
• In a product $AB\mathbf{x}$, the rightmost matrix acts first.

🧪 Practice

1. Verify that $(A + B)(A - B) \ne A^2 - B^2$ for $A = \begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$ and $B = \begin{bmatrix} 1 & 0 \ 1 & 1 \end{bmatrix}$, and expand $(A + B)(A - B)$ correctly.
2. Find all $2 \times 2$ matrices that commute with $\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$. Hint: write the unknown matrix with entries $a, b, c, d$ and compare both products entry by entry.
3. Find a nonzero $2 \times 2$ matrix $A$ with $A^2 = O$. Can such a matrix have a nonzero entry on its diagonal?
4. Interview: A graphics engine stores each object's transform as a product $T R S$ of a translation, rotation, and scale, applied to vertex coordinates as $T R S\mathbf{v}$. A teammate proposes storing $S R T$ instead "because multiplication is associative". What goes wrong, and which property did they confuse? Hint: associativity regroups a product without reordering it; ask what happens to a vertex if it is translated before it is rotated.

Powers of a Matrix

Why it matters. Many processes apply the same linear step over and over: a population model advancing one year at a time, a Markov chain taking one random step, a recurrence such as the Fibonacci numbers producing one term at a time. If one step is $\mathbf{x} \mapsto A\mathbf{x}$, then $k$ steps are $\mathbf{x} \mapsto A^k\mathbf{x}$. Powers of matrices are therefore the algebra of repeated processes, and Chapter 8 shows how eigenvalues predict their long-run behavior.

Definition. For a square matrix $A$ and an integer $k \ge 1$,

$$A^0 = I, \qquad A^k = \underbrace{A A \cdots A}_{k \text{ factors}}.$$

Only square matrices have powers, because $AA$ requires the number of columns of $A$ to equal its number of rows. Associativity makes the grouping irrelevant.

Rules. For nonnegative integers $r$ and $s$:

Rule	Holds?
$A^r A^s = A^{r+s}$	Yes
$(A^r)^s = A^{rs}$	Yes
$A^r A^s = A^s A^r$	Yes: powers of the same matrix commute
$(AB)^k = A^k B^k$	Only if $AB = BA$ in general; $(AB)^2 = ABAB$

Example: Fibonacci numbers. The recurrence $F_{k+1} = F_k + F_{k-1}$ advances the pair $(F_k, F_{k-1})$ by a matrix:

$$\begin{bmatrix} F_{k+1} \ F_k \end{bmatrix} = \begin{bmatrix} 1 & 1 \ 1 & 0 \end{bmatrix} \begin{bmatrix} F_k \ F_{k-1} \end{bmatrix}, \qquad\text{and in fact}\qquad Q^k = \begin{bmatrix} 1 & 1 \ 1 & 0 \end{bmatrix}^k = \begin{bmatrix} F_{k+1} & F_k \ F_k & F_{k-1} \end{bmatrix}.$$

For instance, $Q^{10} = \begin{bmatrix} 89 & 55 \ 55 & 34 \end{bmatrix}$.

Example: counting walks in a graph. The adjacency matrix of a graph with vertices $1, \dots, n$ has $a_{ij} = 1$ if there is an edge from $i$ to $j$, and $0$ otherwise. Then $(A^k){ij}$ is the number of walks of length $k$ from $i$ to $j$. For $k = 2$, $(A^2){ij} = \sum_l a_{il}a_{lj}$ counts the intermediate vertices $l$ with an edge $i \to l$ and an edge $l \to j$. Induction on $k$ gives the general case.

     1 ----- 2              A = [ 0 1 1 0 ]      A^2 = [ 2 1 1 1 ]
      \     /                   [ 1 0 1 0 ]            [ 1 2 1 1 ]
       \   /                    [ 1 1 0 1 ]            [ 1 1 3 0 ]
         3 ----- 4              [ 0 0 1 0 ]            [ 1 1 0 1 ]

   (A^2)_33 = 3: the walks 3-1-3, 3-2-3, 3-4-3
   (A^2)_14 = 1: the walk 1-3-4

Special behavior of powers.

• A matrix with $N^k = O$ for some $k$ is nilpotent, such as $\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$. Numbers cannot do this, apart from $0$.
• A matrix with $P^2 = P$ is idempotent. Projections are idempotent (Section 7.1): projecting twice is the same as projecting once.
• Negative powers $A^{-k} = (A^{-1})^k$ are defined only for invertible matrices (Section 4.3).

Fast exponentiation. Computing $A^k$ by $k - 1$ multiplications is wasteful. Repeated squaring computes $A, A^2, A^4, A^8, \dots$ and multiplies together the powers that correspond to the 1-bits of $k$. For $k = 13 = 1101_2$, $A^{13} = A^8 A^4 A^1$, using 5 multiplications instead of 12. In general, about $2\log_2 k$ multiplications suffice.

import numpy as np

def mat_pow(A, k):
    """A^k for square A and integer k >= 0, by repeated squaring."""
    result = np.eye(A.shape[0], dtype=A.dtype)   # A^0 = I
    base = A.copy()
    while k > 0:
        if k & 1:                                # current lowest bit of k is 1
            result = result @ base
        base = base @ base                       # A, A^2, A^4, A^8, ...
        k >>= 1
    return result

Q = np.array([[1, 1], [1, 0]])
print(mat_pow(Q, 10))                            # [[89 55] [55 34]]
print(np.array_equal(mat_pow(Q, 10), np.linalg.matrix_power(Q, 10)))   # True

# Fixed-size integers overflow silently; use Python integers (dtype=object) for exact large powers
Qbig = np.array([[1, 1], [1, 0]], dtype=object)
print(mat_pow(Qbig, 100)[0, 1])                  # 354224848179261915075 = F_100
print(np.linalg.matrix_power(Q, 100)[0, 1])      # 3736710778780434371: wrong, int64 overflowed

# Walks in the graph 1-2, 1-3, 2-3, 3-4
A = np.array([[0, 1, 1, 0],
              [1, 0, 1, 0],
              [1, 1, 0, 1],
              [0, 0, 1, 0]])
A2 = A @ A
print(A2[2, 2], A2[0, 3])                        # 3 1: walks of length 2 from 3 to 3, and from 1 to 4
print(np.diag(mat_pow(A, 3)) // 2)               # [1 1 1 0]: triangles through each vertex

# A nilpotent matrix: a nonzero matrix whose square is zero
N = np.array([[0, 1], [0, 0]])
print(N @ N)                                     # [[0 0] [0 0]]

A closed walk of length 3 from a vertex goes around a triangle, in one of two directions, so half of $(A^3)_{ii}$ counts the triangles through vertex $i$.

Key takeaways

• Only square matrices have powers; $A^0 = I$ and $A^k$ is the product of $k$ copies of $A$.
• $A^rA^s = A^{r+s}$ and $(A^r)^s = A^{rs}$, but $(AB)^k = A^kB^k$ requires $AB = BA$.
• $A^k\mathbf{x}$ is the result of $k$ steps of the process $\mathbf{x} \mapsto A\mathbf{x}$; $(A^k)_{ij}$ counts walks of length $k$ in a graph.
• Nonzero matrices can be nilpotent ($N^k = O$) or idempotent ($P^2 = P$).
• Repeated squaring computes $A^k$ with about $2\log_2 k$ multiplications.

🧪 Practice

1. Compute $A^2$, $A^3$, and guess a formula for $A^k$, where $A = \begin{bmatrix} 1 & 2 \ 0 & 1 \end{bmatrix}$. Prove the formula by induction.
2. Show that $(AB)^2 = A^2B^2$ fails for $A = \begin{bmatrix} 1 & 1 \ 0 & 0 \end{bmatrix}$ and $B = \begin{bmatrix} 1 & 0 \ 1 & 0 \end{bmatrix}$.
3. For the cycle graph on 4 vertices ($1-2-3-4-1$), write the adjacency matrix $A$ and use $A^3$ to count the walks of length 3 from vertex 1 to vertex 2. List them.
4. Interview: Compute the $10^{18}$-th Fibonacci number modulo $10^9 + 7$. A loop over $k$ is far too slow. What would you do, and how many matrix multiplications does it take? Hint: combine repeated squaring of $Q$ with reducing every entry modulo $10^9 + 7$ after each multiplication.

<a id="42-transpose-and-special-matrices"></a>

4.2 Transpose and Special Matrices

This subchapter introduces the transpose and the families of matrices with special structure (identity, diagonal, triangular, symmetric, Hermitian) that appear constantly in theory and in fast algorithms.

Matrix Transpose and Properties

Why it matters. Many problems need to switch between the row view and the column view of the same data. A data table with one row per customer and one column per feature is naturally read by rows, but statistics such as feature correlations are computed by columns. The transpose exchanges the two views. It also turns the dot product into a matrix product, $\mathbf{u} \cdot \mathbf{v} = \mathbf{u}^T\mathbf{v}$, which is how geometry enters matrix algebra.

Definition. The transpose of an $m \times n$ matrix $A$ is the $n \times m$ matrix $A^T$ with entries

$$(A^T){ij} = a{ji}.$$

Row $i$ of $A$ becomes column $i$ of $A^T$. Geometrically, $A^T$ is $A$ reflected across its main diagonal.

        A (2 x 3)                    A^T (3 x 2)

   [ 1   2   3 ]                     [ 1   4 ]
   [ 4   5   6 ]       ------>       [ 2   5 ]
                                     [ 3   6 ]
   row 1 = (1, 2, 3)                 column 1 = (1, 2, 3)
   a_12 = 2                          (A^T)_21 = 2

The transpose of a column vector is a row vector. This gives a compact way to write vectors in running text, $\mathbf{x} = [,1\ \ 2\ \ 3,]^T$, and to write the dot product as a matrix product:

$$\mathbf{u}^T\mathbf{v} = [,u_1\ \cdots\ u_n,]\begin{bmatrix} v_1 \ \vdots \ v_n \end{bmatrix} = u_1v_1 + \cdots + u_nv_n = \mathbf{u} \cdot \mathbf{v}.$$

The product in the other order, $\mathbf{u}\mathbf{v}^T$, is the $n \times n$ outer product with entries $u_iv_j$.

Properties. For matrices of compatible sizes and a scalar $c$:

Property	Statement
Involution	$(A^T)^T = A$
Sum	$(A + B)^T = A^T + B^T$
Scalar multiple	$(cA)^T = cA^T$
Product (order reverses)	$(AB)^T = B^TA^T$
Several factors	$(A_1A_2\cdots A_k)^T = A_k^T \cdots A_2^TA_1^T$

Why the product rule reverses the order. First, sizes force it. If $A$ is $m \times n$ and $B$ is $n \times p$, then $(AB)^T$ is $p \times m$. The product $A^TB^T$ is $(n \times m)(p \times n)$, which is usually not even defined, while $B^TA^T$ is $(p \times n)(n \times m) = p \times m$. Second, the entries agree:

$$\big((AB)^T\big){ij} = (AB){ji} = \sum_k a_{jk}b_{ki} = \sum_k (B^T){ik}(A^T){kj} = (B^TA^T)_{ij}.$$

The "socks and shoes" analogy applies again: to undo putting on socks and then shoes, you take off the shoes first.

A useful consequence. For any $A$, the dot product of $A\mathbf{x}$ with $\mathbf{y}$ can move $A$ to the other side as a transpose:

$$(A\mathbf{x}) \cdot \mathbf{y} = (A\mathbf{x})^T\mathbf{y} = \mathbf{x}^TA^T\mathbf{y} = \mathbf{x} \cdot (A^T\mathbf{y}).$$

This identity is the reason the transpose appears in least squares (Section 9.5), in the four fundamental subspaces (Section 6.6), and in backpropagation, where gradients flow through a layer $\mathbf{y} = W\mathbf{x}$ by multiplication with $W^T$.

import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])
print(A.T)                                   # [[1 4] [2 5] [3 6]]
print(A.shape, A.T.shape)                    # (2, 3) (3, 2)

rng = np.random.default_rng(1)
B = rng.integers(-5, 5, size=(3, 4))
print(np.array_equal((A @ B).T, B.T @ A.T))  # True: order reverses
# A.T @ B.T would be (3 x 2)(4 x 3): undefined, so NumPy raises a ValueError

# Dot product and outer product via transposes (vectors as n x 1 matrices)
u = np.array([[1], [2], [3]])
v = np.array([[4], [0], [-1]])
print(u.T @ v)                               # [[1]]: 4 + 0 - 3, a 1 x 1 matrix
print(u @ v.T)                               # 3 x 3 outer product, entries u_i * v_j

# (Ax) . y == x . (A^T y)
x = rng.standard_normal(3)
y = rng.standard_normal(2)
print(np.isclose((A @ x) @ y, x @ (A.T @ y)))   # True

# Pitfall: transposing a 1-D NumPy array does nothing
w = np.array([1, 2, 3])
print(w.shape, w.T.shape)                    # (3,) (3,): use w.reshape(-1, 1) for a column

Key takeaways

• $A^T$ swaps rows and columns: $(A^T){ij} = a{ji}$, and an $m \times n$ matrix becomes $n \times m$.
• Transposition respects sums and scalar multiples but reverses products: $(AB)^T = B^TA^T$.
• The dot product is $\mathbf{u}^T\mathbf{v}$ and the outer product is $\mathbf{u}\mathbf{v}^T$.
• $(A\mathbf{x}) \cdot \mathbf{y} = \mathbf{x} \cdot (A^T\mathbf{y})$: the transpose moves a matrix across a dot product.

🧪 Practice

1. For $A = \begin{bmatrix} 2 & -1 \ 0 & 3 \ 1 & 4 \end{bmatrix}$ and $B = \begin{bmatrix} 1 & 2 \ -1 & 0 \end{bmatrix}$, compute $(AB)^T$ and $B^TA^T$ and confirm they are equal.
2. Show that $A^TA$ and $AA^T$ are both defined for any $m \times n$ matrix $A$, give their sizes, and show that the diagonal entries of $A^TA$ are the squared lengths of the columns of $A$.
3. Prove that $(ABC)^T = C^TB^TA^T$ using only the two-factor rule.
4. Interview: In a neural network layer $\mathbf{y} = W\mathbf{x}$ with $W$ of size $256 \times 784$, the gradient of the loss with respect to $\mathbf{x}$ is computed from the gradient $\mathbf{g}$ with respect to $\mathbf{y}$. Why must the formula involve $W^T\mathbf{g}$ rather than $W\mathbf{g}$? Hint: check sizes first, then use $(W\mathbf{x}) \cdot \mathbf{g} = \mathbf{x} \cdot (W^T\mathbf{g})$.

Identity and Zero Matrices

Why it matters. Arithmetic with numbers relies on $1$ (multiplying by it changes nothing) and $0$ (adding it changes nothing). Matrix algebra has the same two anchors: the identity matrix, which is the transformation "do nothing", and the zero matrix. Equations such as $AA^{-1} = I$ and $A + O = A$ cannot even be stated without them.

Identity matrix. The $n \times n$ identity matrix $I_n$ (or simply $I$) has ones on the main diagonal and zeros elsewhere:

$$I_3 = \begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{bmatrix}, \qquad (I_n){ij} = \delta{ij} = \begin{cases} 1 & i = j \ 0 & i \ne j \end{cases}.$$

The symbol $\delta_{ij}$ is the Kronecker delta. The columns of $I_n$ are the standard basis vectors $\mathbf{e}_1, \dots, \mathbf{e}_n$ (Section 2.1).

The identity property. For every $m \times n$ matrix $A$,

$$I_m A = A = A I_n.$$

The two identities have different sizes when $A$ is not square. The proof is one line: $(I_mA){ij} = \sum_k \delta{ik}a_{kj} = a_{ij}$, because only the term $k = i$ survives. As a transformation, $I\mathbf{x} = \mathbf{x}$ for every $\mathbf{x}$.

Extracting columns and rows. Multiplying by a standard basis vector selects part of a matrix:

$$A\mathbf{e}_j = \operatorname{col}_j(A), \qquad \mathbf{e}_i^TA = \operatorname{row}_i(A), \qquad \mathbf{e}_i^TA\mathbf{e}j = a{ij}.$$

This turns statements about individual entries into matrix equations, a trick used in many proofs. For example, if $A\mathbf{x} = B\mathbf{x}$ for every $\mathbf{x}$, then choosing $\mathbf{x} = \mathbf{e}_j$ shows that every column of $A$ equals the corresponding column of $B$, so $A = B$.

Zero matrix. The $m \times n$ zero matrix $O_{m \times n}$ (or $O$, or $0$) has every entry equal to zero. It satisfies

$$A + O = A, \qquad A - A = O, \qquad OA = O, \qquad AO = O,$$

with sizes chosen so that each expression is defined. As a transformation, $O\mathbf{x} = \mathbf{0}$ for every $\mathbf{x}$: it collapses everything to the origin.

Where the analogy with numbers breaks. The number 1 has only two square roots, $\pm 1$. The identity matrix has infinitely many. Every reflection of the plane satisfies $F^2 = I$ because reflecting twice returns every point to where it started:

$$\begin{bmatrix} \cos\theta & \sin\theta \ \sin\theta & -\cos\theta \end{bmatrix}^2 = I \quad\text{for every } \theta.$$

Similarly, as Section 4.1 showed, $O$ has nonzero square roots, such as $\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$.

import numpy as np

I3 = np.eye(3)
A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
print(np.array_equal(I3 @ A, A), np.array_equal(A @ I3, A))   # True True

# Non-square A needs a differently sized identity on each side
M = np.arange(6.).reshape(2, 3)
print(np.array_equal(np.eye(2) @ M, M), np.array_equal(M @ np.eye(3), M))   # True True

# Standard basis vectors select columns, rows, and entries
e = np.eye(3)                              # e[:, j] is e_{j+1}
print(A @ e[:, 1])                         # [ 1. -1.  1.]: column 2 of A
print(e[:, 2] @ A)                         # [-2.  1.  2.]: row 3 of A
print(e[:, 1] @ A @ e[:, 0])               # -3.0: the entry a_21

# Infinitely many square roots of I: every reflection matrix
for theta in [0.3, 1.0, 2.5]:
    F = np.array([[np.cos(theta), np.sin(theta)], [np.sin(theta), -np.cos(theta)]])
    print(np.allclose(F @ F, np.eye(2)))   # True, three times

Z = np.zeros((2, 3))
print(np.array_equal(Z @ np.ones((3, 4)), np.zeros((2, 4))))   # True: OB = O

Key takeaways

• $I_n$ has ones on the diagonal and zeros elsewhere; its columns are $\mathbf{e}_1, \dots, \mathbf{e}n$ and its entries are $\delta{ij}$.
• $I_mA = A = AI_n$ for every $m \times n$ matrix $A$; $I$ is the "do nothing" transformation.
• $A\mathbf{e}_j$ is column $j$ of $A$, $\mathbf{e}_i^TA$ is row $i$, and $\mathbf{e}_i^TA\mathbf{e}j = a{ij}$.
• The zero matrix is the additive identity and annihilates everything it multiplies.
• Unlike the number 1, $I$ has infinitely many square roots, such as all reflections.

🧪 Practice

1. Write $I_4$ and compute $I_4\mathbf{x}$ for $\mathbf{x} = (3, -1, 0, 2)$. Explain why $I_3\mathbf{x}$ is undefined.
2. Use the vectors $\mathbf{e}_j$ to prove: if $A$ is $m \times n$ and $A\mathbf{x} = \mathbf{0}$ for every $\mathbf{x} \in \mathbb{R}^n$, then $A = O$.
3. Find all diagonal $2 \times 2$ matrices $D$ with $D^2 = I$. Then find one $2 \times 2$ matrix with $X^2 = I$ that is not diagonal and not $\pm I$.
4. Interview: Prove that $I$ is the only $n \times n$ matrix $E$ with $EA = A$ for every $n \times n$ matrix $A$. Hint: the property must hold for the particular choice $A = I$.

Diagonal and Scalar Matrices

Why it matters. A diagonal matrix acts on each coordinate separately, scaling coordinate $i$ by $d_i$ and leaving the others alone. This makes diagonal matrices as easy to handle as numbers: their products, powers, and inverses are computed entry by entry. Much of linear algebra, including diagonalization (Chapter 8) and the SVD (Chapter 11), consists of finding coordinates in which a given matrix becomes diagonal.

Definitions. A square matrix $D$ is diagonal if every entry off the main diagonal is zero, $d_{ij} = 0$ for $i \ne j$. It is written

$$D = \operatorname{diag}(d_1, d_2, \dots, d_n) = \begin{bmatrix} d_1 & & 0 \ & \ddots & \ 0 & & d_n \end{bmatrix}.$$

Zeros are allowed on the diagonal itself. A scalar matrix is a diagonal matrix with all diagonal entries equal, $cI = \operatorname{diag}(c, \dots, c)$. The identity and the square zero matrix are both scalar matrices.

Multiplying by a diagonal matrix. The side matters:

Product	Effect	Reason
$D\mathbf{x}$	Scales component $i$ by $d_i$	$(D\mathbf{x})_i = d_ix_i$
$DA$	Scales row $i$ of $A$ by $d_i$	$(DA){ij} = d_ia{ij}$
$AD$	Scales column $j$ of $A$ by $d_j$	$(AD){ij} = a{ij}d_j$
$cIA = AcI$	Scales all of $A$ by $c$	Same as the scalar multiple $cA$

Example. A table $A$ with rows = stores and columns = products holds units sold. Multiplying on the right by $D = \operatorname{diag}(\text{price}_1, \dots, \text{price}_n)$ converts every column from units to revenue. Multiplying on the left by $\operatorname{diag}(1/\text{staff}_1, \dots)$ converts every row to revenue per employee.

Algebra of diagonal matrices. For $D = \operatorname{diag}(d_1, \dots, d_n)$ and $E = \operatorname{diag}(e_1, \dots, e_n)$:

$$DE = ED = \operatorname{diag}(d_1e_1, \dots, d_ne_n), \qquad D^k = \operatorname{diag}(d_1^k, \dots, d_n^k),$$

and $D$ is invertible exactly when every $d_i \ne 0$, with $D^{-1} = \operatorname{diag}(1/d_1, \dots, 1/d_n)$. All diagonal matrices of the same size commute with each other.

Scalar matrices commute with everything. $(cI)A = cA = A(cI)$ for every $A$. The converse is also true: a matrix that commutes with every $n \times n$ matrix must be a scalar matrix (Practice 3).

   D = diag(2, 1/2) acting on the plane:

        y                               y
        |                               |
      1 +-----+                    1/2  +-----------+
        |     |       ------>           |           |
        +-----+---- x                   +-----------+---- x
              1                                     2

   x-coordinates doubled, y-coordinates halved; area scaled by 2 * 1/2 = 1

import numpy as np

d = np.array([2., 3., -1.])
D = np.diag(d)                                  # builds the 3 x 3 diagonal matrix
A = np.arange(1., 10.).reshape(3, 3)            # [[1 2 3] [4 5 6] [7 8 9]]

print(D @ A)                                    # rows scaled: [[2 4 6] [12 15 18] [-7 -8 -9]]
print(A @ D)                                    # columns scaled: [[2 6 -3] [8 15 -6] [14 24 -9]]

# Same results without forming D: n multiplications per row instead of a full matrix product
print(np.array_equal(D @ A, d[:, None] * A))    # True: d[:, None] has shape (3, 1), scales rows
print(np.array_equal(A @ D, A * d))             # True: d broadcasts along rows, scales columns

# Powers and inverses act entry by entry
print(np.allclose(np.linalg.matrix_power(D, 5), np.diag(d ** 5)))   # True
print(np.allclose(np.linalg.inv(D), np.diag(1 / d)))                # True

# Scalar matrices commute with everything; general diagonal matrices do not
print(np.array_equal((4 * np.eye(3)) @ A, A @ (4 * np.eye(3))))    # True
print(np.array_equal(D @ A, A @ D))                                 # False

The broadcasting form in the code matters in practice. Forming $D$ and computing $DA$ costs $n^3$ multiplications, while scaling the rows directly costs $n^2$.

Key takeaways

• A diagonal matrix has zeros off the main diagonal; a scalar matrix is $cI$.
• $DA$ scales the rows of $A$, while $AD$ scales the columns.
• Products, powers, and inverses of diagonal matrices are computed entry by entry; diagonal matrices commute with each other.
• $D$ is invertible exactly when no diagonal entry is zero.
• Scalar matrices are exactly the matrices that commute with every matrix.

🧪 Practice

1. For $D = \operatorname{diag}(3, -2)$ and $A = \begin{bmatrix} 1 & 4 \ 2 & 5 \end{bmatrix}$, compute $DA$, $AD$, $D^4$, and $D^{-1}$.
2. Find a diagonal matrix $D$ such that $DA$ has every row summing to 1, where $A = \begin{bmatrix} 1 & 3 \ 2 & 2 \ 0 & 5 \end{bmatrix}$. (This normalization turns a table of counts into a table of proportions.)
3. Let $A$ be a $2 \times 2$ matrix that commutes with $E_{11} = \begin{bmatrix} 1 & 0 \ 0 & 0 \end{bmatrix}$ and $E_{12} = \begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$. Show that $A$ is a scalar matrix.
4. Interview: A feature matrix $X$ has shape $(10^6, 50)$, and you need to standardize each column by dividing by its standard deviation. A colleague writes X @ np.diag(1 / std). Is the result correct? Is it efficient? Hint: compare the cost of the matrix product with that of broadcasting X / std, and consider memory as well as time.

Triangular Matrices

Why it matters. Triangular systems are the easy systems: they are solved one unknown at a time by substitution, in about $n^2$ operations instead of the $n^3$ needed for a general system. Gaussian elimination is a method for reaching triangular form, and the LU factorization (Section 4.4) writes any matrix in terms of triangular ones. Recognizing and preserving triangular structure is one of the main sources of speed in numerical linear algebra.

Definitions. A square matrix $A$ is

• upper triangular if $a_{ij} = 0$ whenever $i > j$ (all entries below the diagonal are zero);
• lower triangular if $a_{ij} = 0$ whenever $i < j$ (all entries above the diagonal are zero);
• unit (upper or lower) triangular if, in addition, every diagonal entry is $1$;
• strictly triangular if the diagonal is zero as well.

A matrix is diagonal exactly when it is both upper and lower triangular.

     upper triangular         lower triangular        unit lower triangular

     [ *  *  *  * ]           [ *  0  0  0 ]          [ 1  0  0  0 ]
     [ 0  *  *  * ]           [ *  *  0  0 ]          [ *  1  0  0 ]
     [ 0  0  *  * ]           [ *  *  *  0 ]          [ *  *  1  0 ]
     [ 0  0  0  * ]           [ *  *  *  * ]          [ *  *  *  1 ]

     (* = any value, possibly zero)

Closure properties. For $n \times n$ upper triangular $U_1$ and $U_2$ (and likewise for lower triangular):

Operation	Result
$U_1 + U_2$, $cU_1$	Upper triangular
$U_1U_2$	Upper triangular, with diagonal entries $(U_1){ii}(U_2){ii}$
$U_1^{-1}$ (if it exists)	Upper triangular, with diagonal entries $1/(U_1)_{ii}$
$U_1^T$	Lower triangular
Product of unit triangular matrices	Unit triangular

Why the product is upper triangular. For $i > j$, $(U_1U_2){ij} = \sum_k (U_1){ik}(U_2){kj}$. A term can be nonzero only if $k \ge i$ (otherwise $(U_1){ik} = 0$) and $k \le j$ (otherwise $(U_2)_{kj} = 0$). With $i > j$, no $k$ satisfies both, so every term is zero. For $i = j$, only $k = i$ survives, giving the product of diagonal entries.

Invertibility. A triangular matrix is invertible exactly when every diagonal entry is nonzero. The diagonal entries are the pivots: if they are all nonzero, back substitution never divides by zero and $U\mathbf{x} = \mathbf{b}$ has exactly one solution for every $\mathbf{b}$. If some diagonal entry is zero, there are fewer than $n$ pivots and the system cannot have a unique solution (Section 4.3).

Forward substitution. A lower triangular system $L\mathbf{x} = \mathbf{b}$ is solved from the top down, mirroring back substitution:

$$x_i = \frac{1}{l_{ii}}\Big(b_i - \sum_{j < i} l_{ij}x_j\Big), \qquad i = 1, 2, \dots, n.$$

Step $i$ uses $i - 1$ multiplications, so the total is about $n^2/2$ multiplications.

Example.

$$\begin{bmatrix} 2 & 0 & 0 \ 1 & 3 & 0 \ -1 & 2 & 4 \end{bmatrix}\begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix} = \begin{bmatrix} 4 \ 5 \ 6 \end{bmatrix}: \qquad x_1 = 2,\quad x_2 = \frac{5 - 2}{3} = 1,\quad x_3 = \frac{6 + 2 - 2}{4} = \frac{3}{2}.$$

import numpy as np
from scipy.linalg import solve_triangular

def forward_substitution(L, b):
    """Solve Lx = b for lower triangular L with nonzero diagonal, in O(n^2) operations."""
    n = len(b)
    x = np.zeros(n)
    for i in range(n):
        x[i] = (b[i] - L[i, :i] @ x[:i]) / L[i, i]   # only x_1 .. x_{i-1} are needed, all already known
    return x

L = np.array([[2., 0, 0], [1, 3, 0], [-1, 2, 4]])
b = np.array([4., 5, 6])
print(forward_substitution(L, b))              # [2.  1.  1.5]
print(solve_triangular(L, b, lower=True))      # [2.  1.  1.5]: the library routine

# Products and inverses stay triangular; diagonals multiply
U1 = np.triu(np.arange(1., 10.).reshape(3, 3))  # np.triu zeroes everything below the diagonal
U2 = np.triu(np.ones((3, 3))) * 2
P = U1 @ U2
print(np.allclose(P, np.triu(P)))              # True: the product is upper triangular
print(np.diag(P), np.diag(U1) * np.diag(U2))   # [ 2. 10. 18.] [ 2. 10. 18.]
print(np.allclose(np.linalg.inv(U1), np.triu(np.linalg.inv(U1))))   # True

Key takeaways

• Upper triangular matrices have zeros below the diagonal, lower triangular matrices have zeros above it, and unit triangular matrices also have ones on the diagonal.
• Sums, products, and inverses of upper (lower) triangular matrices are upper (lower) triangular; the transpose switches the type.
• The diagonal of a product of triangular matrices is the product of their diagonals.
• A triangular matrix is invertible exactly when all diagonal entries are nonzero.
• Triangular systems are solved by forward or back substitution in about $n^2$ operations instead of $n^3$.

🧪 Practice

1. Solve $\begin{bmatrix} 1 & 0 & 0 \ 2 & 1 & 0 \ -3 & 4 & 1 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 1 \ 0 \ 2 \end{bmatrix}$ by forward substitution.
2. Compute the product of $\begin{bmatrix} 1 & 2 \ 0 & 3 \end{bmatrix}$ and $\begin{bmatrix} 4 & -1 \ 0 & 5 \end{bmatrix}$ in both orders and confirm that both are upper triangular with the same diagonal. Are the products equal?
3. Show that a strictly upper triangular $3 \times 3$ matrix $N$ satisfies $N^3 = O$. Hint: track how the band of zeros grows with each multiplication.
4. Interview: You must solve $T\mathbf{x} = \mathbf{b}$ where $T$ is a $5000 \times 5000$ upper triangular matrix. A colleague calls a general solver, which spends most of its time factoring $T$. Estimate the speedup from calling a triangular solver instead. Hint: compare about $\tfrac23 n^3$ operations with about $n^2$ for $n = 5000$.

Symmetric and Skew-Symmetric Matrices

Why it matters. Many matrices describe a relationship between pairs that does not depend on the order of the pair: the distance between cities $i$ and $j$, the covariance between features $i$ and $j$, the stiffness coupling between joints $i$ and $j$, the second partial derivatives $\partial^2 f / \partial x_i \partial x_j$. Such matrices equal their own transposes. Symmetric matrices have the best-behaved theory in linear algebra (real eigenvalues and perpendicular eigenvectors, Chapter 10), so recognizing symmetry pays off.

Definitions. A square matrix $A$ is

• symmetric if $A^T = A$, that is, $a_{ij} = a_{ji}$ for all $i, j$;
• skew-symmetric (or antisymmetric) if $A^T = -A$, that is, $a_{ij} = -a_{ji}$ for all $i, j$.

Setting $i = j$ in the skew-symmetric condition gives $a_{ii} = -a_{ii}$, so the diagonal of a skew-symmetric matrix is zero.

       symmetric                     skew-symmetric

    [  4   1  -2 ]                  [  0   3  -1 ]
    [  1   0   5 ]                  [ -3   0   2 ]
    [ -2   5   7 ]                  [  1  -2   0 ]

    mirror image across the          mirror image with sign flipped;
    diagonal                         diagonal forced to 0

Sources of symmetric matrices.

• For any $m \times n$ matrix $A$, the products $A^TA$ ($n \times n$) and $AA^T$ ($m \times m$) are symmetric: $(A^TA)^T = A^T(A^T)^T = A^TA$. The entry $(A^TA)_{ij}$ is the dot product of columns $i$ and $j$ of $A$, so $A^TA$ is the Gram matrix of the columns.
• For a square $A$, $A + A^T$ is symmetric and $A - A^T$ is skew-symmetric.

Symmetric-skew decomposition. Every square matrix splits uniquely into a symmetric part and a skew-symmetric part:

$$A = \underbrace{\tfrac12(A + A^T)}{\text{symmetric}} + \underbrace{\tfrac12(A - A^T)}{\text{skew-symmetric}}.$$

This mirrors the split of a function into even and odd parts. For uniqueness, if $A = S + K$ with $S$ symmetric and $K$ skew-symmetric, then $A^T = S - K$, and adding and subtracting the two equations recovers the formulas above.

Products need care. Sums and scalar multiples of symmetric matrices are symmetric, but products usually are not. If $A$ and $B$ are symmetric, $(AB)^T = B^TA^T = BA$, so $AB$ is symmetric exactly when $AB = BA$.

Skew-symmetric matrices in geometry. The cross product with a fixed vector $\mathbf{a} \in \mathbb{R}^3$ (Section 2.4) is multiplication by a skew-symmetric matrix:

$$\mathbf{a} \times \mathbf{b} = [\mathbf{a}]\times,\mathbf{b}, \qquad [\mathbf{a}]\times = \begin{bmatrix} 0 & -a_3 & a_2 \ a_3 & 0 & -a_1 \ -a_2 & a_1 & 0 \end{bmatrix}.$$

For any skew-symmetric $K$, $\mathbf{x}^TK\mathbf{x} = 0$ for every $\mathbf{x}$, because this $1 \times 1$ number equals its own transpose $\mathbf{x}^TK^T\mathbf{x} = -\mathbf{x}^TK\mathbf{x}$. For the cross product, this says $\mathbf{x} \cdot (\mathbf{a} \times \mathbf{x}) = 0$: the cross product is perpendicular to its factors. Skew-symmetric matrices are also the infinitesimal rotations of Section 9.4.

import numpy as np

A = np.array([[1., 2, 0], [4, 3, -1], [2, 5, 6]])
S = (A + A.T) / 2                              # symmetric part
K = (A - A.T) / 2                              # skew-symmetric part
print(np.array_equal(S, S.T), np.array_equal(K, -K.T), np.array_equal(S + K, A))   # True True True
print(np.diag(K))                              # [0. 0. 0.]

# A^T A is symmetric for any rectangular A: it is the Gram matrix of the columns
X = np.random.default_rng(2).standard_normal((5, 3))
G = X.T @ X
print(np.allclose(G, G.T), np.isclose(G[0, 1], X[:, 0] @ X[:, 1]))   # True True

# Products of symmetric matrices are usually not symmetric
S2 = np.array([[1., 0, 0], [0, 2, 0], [0, 0, 3]])
print(np.allclose(S @ S2, (S @ S2).T))         # False

# The cross product as a skew-symmetric matrix
def cross_matrix(a):
    return np.array([[0, -a[2], a[1]],
                     [a[2], 0, -a[0]],
                     [-a[1], a[0], 0]])

a, b = np.array([1., 2, 3]), np.array([-2., 0, 5])
print(cross_matrix(a) @ b, np.cross(a, b))     # [10. -11.  4.] [10. -11.  4.]
print(b @ cross_matrix(a) @ b)                 # 0.0: x^T K x = 0 for skew-symmetric K

Key takeaways

• Symmetric means $A^T = A$; skew-symmetric means $A^T = -A$, which forces a zero diagonal.
• $A^TA$ and $AA^T$ are symmetric for every matrix $A$; $A^TA$ holds the dot products of the columns of $A$.
• Every square matrix is uniquely the sum of a symmetric part $\tfrac12(A + A^T)$ and a skew-symmetric part $\tfrac12(A - A^T)$.
• The product of two symmetric matrices is symmetric exactly when they commute.
• $\mathbf{x}^TK\mathbf{x} = 0$ for skew-symmetric $K$; the cross product $\mathbf{a} \times \mathbf{b}$ is $[\mathbf{a}]_\times\mathbf{b}$.

🧪 Practice

1. Decompose $A = \begin{bmatrix} 2 & 5 \ -1 & 3 \end{bmatrix}$ into its symmetric and skew-symmetric parts.
2. Prove that if $A$ is symmetric and $B$ is any matrix of compatible size, then $B^TAB$ is symmetric.
3. Show that the only matrix that is both symmetric and skew-symmetric is the zero matrix. How many free parameters does a general $n \times n$ symmetric matrix have? A skew-symmetric one?
4. Interview: A machine learning model computes a Hessian matrix in floating point, and the result is not exactly symmetric because of rounding. Downstream code requires symmetry. What is the simplest principled fix, and why is it the closest symmetric matrix? Hint: use the symmetric-skew decomposition, and note that the symmetric and skew-symmetric parts are orthogonal under the inner product $\langle A, B\rangle = \sum_{ij} a_{ij}b_{ij}$.

Conjugate Transpose and Hermitian Matrices

Why it matters. Complex matrices are essential in signal processing (the Fourier transform), quantum mechanics (states and observables), and in eigenvalue problems, where a real matrix can have complex eigenvalues (Section 8.1). Over $\mathbb{C}$, the plain transpose gives the wrong notion of length: for $\mathbf{z} = (1, i)$, $\mathbf{z}^T\mathbf{z} = 1 + i^2 = 0$, even though $\mathbf{z} \ne \mathbf{0}$. Conjugating the entries fixes this, and the resulting conjugate transpose takes over the role of the transpose for complex matrices.

Definition. The conjugate transpose (or Hermitian adjoint) of a complex $m \times n$ matrix $A$ is the $n \times m$ matrix

$$A^* = \overline{A}^{,T}, \qquad (A^*){ij} = \overline{a{ji}}.$$

It is also written $A^H$ (common in engineering) or $A^\dagger$ (common in physics). Conjugation and transposition commute, so the order of the two steps does not matter. For a real matrix, $A^* = A^T$.

Example.

$$A = \begin{bmatrix} 1 + i & 2 \ 3i & 4 - 2i \end{bmatrix} \qquad\Longrightarrow\qquad A^* = \begin{bmatrix} 1 - i & -3i \ 2 & 4 + 2i \end{bmatrix}.$$

Length and inner product. For $\mathbf{z} \in \mathbb{C}^n$,

$$\mathbf{z}^*\mathbf{z} = \overline{z_1}z_1 + \cdots + \overline{z_n}z_n = |z_1|^2 + \cdots + |z_n|^2 \ge 0,$$

with equality only for $\mathbf{z} = \mathbf{0}$ (Section 1.2). So $|\mathbf{z}| = \sqrt{\mathbf{z}^\mathbf{z}}$ is a genuine length, and $\mathbf{w}^\mathbf{z}$ is the standard complex inner product (Section 9.1). For $\mathbf{z} = (1, i)$, $\mathbf{z}^*\mathbf{z} = 1 + 1 = 2$.

Properties. They parallel the transpose rules, except that scalars are conjugated:

Property	Statement
Involution	$(A^)^ = A$
Sum	$(A + B)^* = A^* + B^*$
Scalar multiple	$(cA)^* = \overline{c},A^*$
Product	$(AB)^* = B^A^$

Hermitian and skew-Hermitian matrices. A square complex matrix is

• Hermitian if $A^* = A$, that is, $a_{ij} = \overline{a_{ji}}$;
• skew-Hermitian if $A^* = -A$.

The diagonal entries of a Hermitian matrix satisfy $a_{ii} = \overline{a_{ii}}$, so they are real. A real symmetric matrix is Hermitian, so Hermitian matrices are the complex generalization of symmetric ones. Their key property: for Hermitian $A$, the number $\mathbf{z}^*A\mathbf{z}$ is real for every $\mathbf{z}$. Indeed, this $1 \times 1$ matrix equals its own conjugate transpose, $(\mathbf{z}^A\mathbf{z})^ = \mathbf{z}^A^\mathbf{z} = \mathbf{z}^*A\mathbf{z}$, and a complex number equal to its conjugate is real. This is why quantum observables, which must have real measured values, are modeled by Hermitian matrices.

Real matrices	Complex analogue
Transpose $A^T$	Conjugate transpose $A^*$
Symmetric, $A^T = A$	Hermitian, $A^* = A$
Skew-symmetric, $A^T = -A$	Skew-Hermitian, $A^* = -A$
Orthogonal, $Q^TQ = I$ (Section 9.4)	Unitary, $U^*U = I$ (Section 9.4)

Example. The Pauli matrices of quantum mechanics are Hermitian:

$$\sigma_x = \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}, \qquad \sigma_y = \begin{bmatrix} 0 & -i \ i & 0 \end{bmatrix}, \qquad \sigma_z = \begin{bmatrix} 1 & 0 \ 0 & -1 \end{bmatrix}.$$

They also illustrate non-commutativity: $\sigma_x\sigma_y = i\sigma_z$, while $\sigma_y\sigma_x = -i\sigma_z$.

import numpy as np

A = np.array([[1 + 1j, 2], [3j, 4 - 2j]])
print(A.conj().T)                     # [[1.-1.j 0.-3.j] [2.-0.j 4.+2.j]]

# Pitfall: .T alone does not conjugate, so z^T z can vanish for a nonzero complex vector
z = np.array([1, 1j])
print(z @ z)                          # 0j: wrong as a squared length
print(z.conj() @ z)                   # (2+0j): the correct squared length
print(np.vdot(z, z))                  # (2+0j): vdot conjugates its first argument

# Pauli matrices: Hermitian, with real diagonals, and they do not commute
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]])
sz = np.array([[1, 0], [0, -1]], dtype=complex)
print(all(np.allclose(s, s.conj().T) for s in (sx, sy, sz)))   # True
print(np.allclose(sx @ sy, 1j * sz), np.allclose(sy @ sx, -1j * sz))   # True True

# z^* H z is real for Hermitian H, for any complex vector z
rng = np.random.default_rng(3)
M = rng.standard_normal((3, 3)) + 1j * rng.standard_normal((3, 3))
H = M + M.conj().T                    # M + M^* is always Hermitian
w = rng.standard_normal(3) + 1j * rng.standard_normal(3)
print(abs((w.conj() @ H @ w).imag) < 1e-12)   # True: the imaginary part is only rounding error

Key takeaways

• The conjugate transpose $A^* = \overline{A}^{,T}$ replaces the transpose for complex matrices; for real matrices they coincide.
• $\mathbf{z}^*\mathbf{z} = \sum |z_i|^2$ is a genuine squared length, while $\mathbf{z}^T\mathbf{z}$ can be zero or complex.
• $(AB)^* = B^A^$ and $(cA)^* = \overline{c}A^*$.
• Hermitian matrices ($A^* = A$) have real diagonals, and $\mathbf{z}^*A\mathbf{z}$ is real for every $\mathbf{z}$.
• In NumPy, use A.conj().T and np.vdot; A.T and @ alone do not conjugate.

🧪 Practice

1. Compute $A^*$ and $A^*A$ for $A = \begin{bmatrix} 1 & i \ 2 - i & 3 \end{bmatrix}$, and check that $A^*A$ is Hermitian.
2. Which of these are Hermitian, skew-Hermitian, or neither? $\begin{bmatrix} 2 & 1 - i \ 1 + i & -5 \end{bmatrix}$, $\ \begin{bmatrix} i & 2 \ -2 & 3i \end{bmatrix}$, $\ \begin{bmatrix} 1 & i \ i & 1 \end{bmatrix}$.
3. Show that $A$ is Hermitian if and only if $iA$ is skew-Hermitian, and conclude that the diagonal entries of a skew-Hermitian matrix are purely imaginary.
4. Interview: A signal-processing routine computes the energy of a complex signal $\mathbf{z}$ as np.dot(z, z) and occasionally returns negative or complex values. What is the bug, and what is the fix? Hint: compare $\sum z_k^2$ with $\sum |z_k|^2$ for $z_k = i$.

<a id="43-matrix-inverses"></a>

4.3 Matrix Inverses

This subchapter asks when a matrix can be "divided by", that is, when its action can be undone, and develops the formulas, algorithms, and equivalent conditions for invertibility.

Definition of the Inverse

Why it matters. The scalar equation $ax = b$ is solved by multiplying by $a^{-1}$, provided $a \ne 0$. The matrix equation $A\mathbf{x} = \mathbf{b}$ invites the same move. Viewing $A$ as a transformation, the inverse is the transformation that undoes it: if $A$ encodes a message, rotates an image, or converts between coordinate systems, $A^{-1}$ decodes, rotates back, or converts back. Unlike nonzero numbers, many nonzero matrices have no inverse, and telling the two kinds apart is one of the central problems of the subject.

Definition. An $n \times n$ matrix $A$ is invertible (or nonsingular) if there is an $n \times n$ matrix $B$ with

$$AB = I_n \qquad\text{and}\qquad BA = I_n.$$

Such a $B$ is called an inverse of $A$. A square matrix with no inverse is singular.

The inverse is unique. Suppose $B$ and $C$ are both inverses of $A$. Then, using associativity,

$$B = BI = B(AC) = (BA)C = IC = C.$$

So we can speak of the inverse and write it $A^{-1}$. The notation $\tfrac{1}{A}$ or $\tfrac{B}{A}$ is never used, because $BA^{-1}$ and $A^{-1}B$ are usually different, and a fraction would not say which one is meant.

Example. For $A = \begin{bmatrix} 2 & 1 \ 5 & 3 \end{bmatrix}$ and $B = \begin{bmatrix} 3 & -1 \ -5 & 2 \end{bmatrix}$:

$$AB = \begin{bmatrix} 6 - 5 & -2 + 2 \ 15 - 15 & -5 + 6 \end{bmatrix} = I, \qquad BA = \begin{bmatrix} 6 - 5 & 3 - 3 \ -10 + 10 & -5 + 6 \end{bmatrix} = I,$$

so $B = A^{-1}$.

Nonzero does not mean invertible. The matrix $S = \begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}$ has no inverse. For any $B$, both rows of $SB$ equal $\operatorname{row}_1(B) + \operatorname{row}_2(B)$, so $SB$ has two equal rows and can never be $I$. Geometrically, $S$ sends every vector $(x, y)$ to $(x + y, x + y)$, which lies on the line $y = x$. Once the plane has been flattened onto a line, there is no way to recover which point a given output came from.

   invertible A: reshapes the plane, reversibly      singular S: flattens the plane onto a line

     +----+        A        /\                        +----+        S          /
     |    |   ---------->  /  \                       |    |   ---------->     /   each output point
     +----+   <----------  \  /                       +----+                 /    comes from many inputs
                 A^-1       \/                                              /     (no way back)

Only square matrices have two-sided inverses. A non-square matrix may have a left inverse ($CA = I$) or a right inverse ($AD = I$), but never both. For example, a $3 \times 2$ matrix maps $\mathbb{R}^2$ into $\mathbb{R}^3$, and its outputs fill at most a plane, so it cannot reach every vector of $\mathbb{R}^3$, which an inverse would require. One-sided inverses appear again with least squares and the pseudoinverse (Sections 9.5 and 11.2).

One check suffices for square matrices. For square $A$ and $B$, $AB = I$ already implies $BA = I$. This is not obvious, and the proof uses the Invertible Matrix Theorem at the end of this subchapter. In practice, it means verifying one product is enough.

import numpy as np

A = np.array([[2., 1], [5, 3]])
B = np.array([[3., -1], [-5, 2]])
I = np.eye(2)
print(np.array_equal(A @ B, I), np.array_equal(B @ A, I))   # True True: B is the inverse of A
print(np.linalg.inv(A))                                     # [[ 3. -1.] [-5.  2.]]

# A nonzero singular matrix: S B always has two equal rows, so it can never be I
S = np.array([[1., 1], [1, 1]])
rng = np.random.default_rng(4)
Bt = rng.standard_normal((2, 2))
print(np.allclose((S @ Bt)[0], (S @ Bt)[1]))                # True
try:
    np.linalg.inv(S)
except np.linalg.LinAlgError as e:
    print("LinAlgError:", e)                                # LinAlgError: Singular matrix

# A^-1 undoes A
x = np.array([4., -1])
print(np.linalg.inv(A) @ (A @ x))                           # [ 4. -1.]

Key takeaways

• $A^{-1}$ is the matrix with $AA^{-1} = A^{-1}A = I$; only square matrices can have one.
• If an inverse exists, it is unique.
• Many nonzero matrices are singular: they collapse space in a way that cannot be undone.
• Never write $B/A$; write $BA^{-1}$ or $A^{-1}B$, which are usually different.
• For square matrices, $AB = I$ implies $BA = I$, so checking one product is enough.

🧪 Practice

1. Verify that $\begin{bmatrix} 1 & 0 & 2 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{bmatrix}$ and $\begin{bmatrix} 1 & 0 & -2 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{bmatrix}$ are inverses of each other, and describe in words what each does to a vector.
2. Show that $\begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}$ is singular by finding a nonzero vector $\mathbf{x}$ with $A\mathbf{x} = \mathbf{0}$, and explain why such an $\mathbf{x}$ rules out an inverse.
3. Find a left inverse of $A = \begin{bmatrix} 1 & 0 \ 0 & 1 \ 0 & 0 \end{bmatrix}$, i.e. a $2 \times 3$ matrix $C$ with $CA = I_2$. Is it unique? Show that $AC \ne I_3$.
4. Interview: Prove that if a square matrix $A$ satisfies $A^2 - 3A + 2I = O$, then $A$ is invertible, and express $A^{-1}$ as a polynomial in $A$. Hint: rearrange to $A \cdot (\text{something}) = I$.

Invertibility of 2x2 Matrices

Why it matters. For $2 \times 2$ matrices, invertibility is decided by a single number, and the inverse has a closed-form formula. This case is worth memorizing: it appears constantly in hand calculations, in geometry, and in proofs, and it previews the determinant, which plays the same role in every dimension (Chapter 5).

Theorem. Let $A = \begin{bmatrix} a & b \ c & d \end{bmatrix}$. Then $A$ is invertible if and only if $ad - bc \ne 0$, in which case

$$A^{-1} = \frac{1}{ad - bc}\begin{bmatrix} d & -b \ -c & a \end{bmatrix}.$$

The number $\det A = ad - bc$ is the determinant of $A$. The recipe: swap the diagonal entries, negate the off-diagonal entries, and divide by the determinant.

Proof. Let $\operatorname{adj} A = \begin{bmatrix} d & -b \ -c & a \end{bmatrix}$. Direct multiplication gives

$$A \operatorname{adj} A = \begin{bmatrix} ad - bc & -ab + ba \ cd - dc & -cb + da \end{bmatrix} = (ad - bc),I,$$

and similarly $(\operatorname{adj} A)A = (ad - bc)I$. If $ad - bc \ne 0$, dividing by it shows that $\frac{1}{ad - bc}\operatorname{adj} A$ is the inverse. If $ad - bc = 0$, then $A \operatorname{adj} A = O$. When $A \ne O$, the matrix $\operatorname{adj} A$ is also nonzero and has a nonzero column $\mathbf{v}$ with $A\mathbf{v} = \mathbf{0}$. An invertible $A$ would give $\mathbf{v} = A^{-1}A\mathbf{v} = \mathbf{0}$, a contradiction. The zero matrix is singular as well. $\blacksquare$

Example. For $A = \begin{bmatrix} 4 & 7 \ 2 & 6 \end{bmatrix}$: $\det A = 24 - 14 = 10$, so

$$A^{-1} = \frac{1}{10}\begin{bmatrix} 6 & -7 \ -2 & 4 \end{bmatrix} = \begin{bmatrix} 0.6 & -0.7 \ -0.2 & 0.4 \end{bmatrix}.$$

Geometric meaning. The columns $(a, c)$ and $(b, d)$ of $A$ are the images of $\mathbf{e}_1$ and $\mathbf{e}_2$. They span a parallelogram whose signed area is $ad - bc$ (Section 2.4). The determinant is zero exactly when the columns are parallel. Then $A$ squashes the whole plane onto a line (or a point), which cannot be undone.

     det != 0                                det = 0

          (b,d)                                    (b,d)
           ^ ______                                 ^
           |/     /     the unit square              /   columns on one line:
           /     /      maps to a genuine           /    the unit square maps to a
          /_____/       parallelogram              /     flat segment of area 0
          +------> (a,c)                          +------> (a,c)

Matrix	$ad - bc$	Invertible?	Geometry
$\begin{bmatrix} 2 & 0 \ 0 & 3 \end{bmatrix}$	6	Yes	Stretch by 2 and 3
$\begin{bmatrix} 0 & -1 \ 1 & 0 \end{bmatrix}$	1	Yes	Rotation by $90^\circ$, inverse rotates back
$\begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}$	0	No	Columns parallel: plane collapses onto a line
$\begin{bmatrix} 1 & k \ 0 & 1 \end{bmatrix}$	1	Yes	Shear; inverse is the shear by $-k$

from fractions import Fraction as Fr

def inverse_2x2(A):
    """Exact inverse of a 2 x 2 matrix [[a, b], [c, d]], or None if it is singular."""
    (a, b), (c, d) = A
    det = Fr(a) * d - Fr(b) * c
    if det == 0:
        return None                                   # columns are parallel: no inverse
    return [[ d / det, -b / det],                     # swap the diagonal, negate the off-diagonal,
            [-c / det,  a / det]]                     # divide everything by the determinant

def show(M):
    return None if M is None else [[str(x) for x in row] for row in M]

print(show(inverse_2x2([[4, 7], [2, 6]])))   # [['3/5', '-7/10'], ['-1/5', '2/5']]
print(show(inverse_2x2([[2, 1], [5, 3]])))   # [['3', '-1'], ['-5', '2']]
print(show(inverse_2x2([[1, 2], [2, 4]])))   # None

# A family of matrices: invertible except at the values of k where the determinant vanishes
for k in [-2, 0, 2, 3]:
    M = [[k, 2], [2, k]]                     # det = k^2 - 4
    print(k, inverse_2x2(M) is not None)     # -2 False, 0 True, 2 False, 3 True

Key takeaways

• $\begin{bmatrix} a & b \ c & d \end{bmatrix}$ is invertible if and only if $ad - bc \ne 0$.
• The inverse is $\frac{1}{ad - bc}\begin{bmatrix} d & -b \ -c & a \end{bmatrix}$: swap, negate, divide.
• $ad - bc$ is the signed area of the parallelogram spanned by the columns; it is zero exactly when the columns are parallel.
• Always check the answer by multiplying: $AA^{-1}$ should be $I$.

🧪 Practice

1. Find the inverses of $\begin{bmatrix} 3 & 5 \ 1 & 2 \end{bmatrix}$ and $\begin{bmatrix} 2 & -4 \ 1 & 3 \end{bmatrix}$, and verify each by multiplication.
2. For which values of $t$ is $\begin{bmatrix} t & 3 \ 4 & t + 1 \end{bmatrix}$ singular?
3. Use the formula to show that the inverse of the rotation matrix $R_\theta = \begin{bmatrix} \cos\theta & -\sin\theta \ \sin\theta & \cos\theta \end{bmatrix}$ is $R_{-\theta}$, and explain why this is geometrically obvious.
4. Interview: A colleague tests whether a matrix is close to singular by checking whether $|\det A| < 10^{-10}$. Show that this test is misleading by evaluating it on $A = 10^{-6}I_2$ and on $A = \begin{bmatrix} 1 & 1 \ 1 & 1 + 10^{-9} \end{bmatrix} \cdot 10^{6}$. What should be used instead? Hint: the determinant changes when you rescale the matrix, but being "nearly singular" should not; see the condition number (Section 11.3).

Computing Inverses by Row Reduction

Why it matters. The $2 \times 2$ formula does not extend conveniently to larger matrices. Row reduction gives a single procedure that, for any square matrix, either computes the inverse or proves that none exists. It also explains what an inverse is: the collection of solutions of $n$ linear systems.

The idea. Finding $X = A^{-1}$ means solving $AX = I$. By the column view of matrix multiplication, this is $n$ separate systems, one for each column of $X$:

$$A\mathbf{x}_1 = \mathbf{e}_1, \quad A\mathbf{x}_2 = \mathbf{e}_2, \quad \dots, \quad A\mathbf{x}_n = \mathbf{e}_n.$$

All of them have the same coefficient matrix, so they can be solved together by augmenting with all the right-hand sides at once (Section 3.2), which is the augmented matrix $[,A \mid I,]$.

Algorithm.

1. Form the $n \times 2n$ matrix $[,A \mid I,]$.
2. Row reduce it to reduced row echelon form (Gauss-Jordan elimination, Section 3.3).
3. If the left block becomes $I$, the right block is $A^{-1}$: $\ [,A \mid I,] \sim [,I \mid A^{-1},]$.
4. If the left block has a row of zeros (fewer than $n$ pivots), $A$ is singular. Stop.

Example. Invert $K = \begin{bmatrix} 1 & 2 & 3 \ 0 & 1 & 4 \ 5 & 6 & 0 \end{bmatrix}$:

  [ 1  2   3 |  1  0  0 ]                            [ 1  2   3 |  1  0  0 ]
  [ 0  1   4 |  0  1  0 ]   R3 <- R3 - 5 R1  --->    [ 0  1   4 |  0  1  0 ]
  [ 5  6   0 |  0  0  1 ]                            [ 0 -4 -15 | -5  0  1 ]

                            R3 <- R3 + 4 R2  --->    [ 1  2   3 |  1  0  0 ]
                                                     [ 0  1   4 |  0  1  0 ]
                                                     [ 0  0   1 | -5  4  1 ]   <- echelon form, 3 pivots

  R2 <- R2 - 4 R3,  R1 <- R1 - 3 R3   --->           [ 1  2   0 |  16 -12 -3 ]
                                                     [ 0  1   0 |  20 -15 -4 ]
                                                     [ 0  0   1 |  -5   4  1 ]

  R1 <- R1 - 2 R2                     --->           [ 1  0   0 | -24  18  5 ]
                                                     [ 0  1   0 |  20 -15 -4 ]
                                                     [ 0  0   1 |  -5   4  1 ]
                                                          I          K^-1

So $K^{-1} = \begin{bmatrix} -24 & 18 & 5 \ 20 & -15 & -4 \ -5 & 4 & 1 \end{bmatrix}$. Checking one entry of $KK^{-1}$: $\operatorname{row}_1(K) \cdot \operatorname{col}_1(K^{-1}) = -24 + 40 - 15 = 1$.

Singular example. For $\begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}$, the step $R_2 \leftarrow R_2 - 2R_1$ turns $[,A \mid I,]$ into $\left[\begin{array}{cc|cc} 1 & 2 & 1 & 0 \ 0 & 0 & -2 & 1 \end{array}\right]$. The zero row on the left means only one pivot, so the matrix is singular.

Why the right block is $A^{-1}$. Each row operation can be performed by multiplying on the left by an invertible matrix (Section 4.4). If the operations $E_1, \dots, E_k$ turn $A$ into $I$, then $E_k \cdots E_1 A = I$, so $E_k \cdots E_1 = A^{-1}$. The same operations applied to the right block turn $I$ into $E_k \cdots E_1 I = A^{-1}$.

Cost. Inverting an $n \times n$ matrix costs about $2n^3$ arithmetic operations, roughly three times as much as one elimination for a single system (Section 3.3).

from fractions import Fraction as Fr

def inverse_gauss_jordan(A):
    """Invert a square matrix exactly by row reducing [A | I]; return None if A is singular."""
    n = len(A)
    M = [[Fr(x) for x in row] + [Fr(int(i == j)) for j in range(n)]   # build [A | I]
         for i, row in enumerate(A)]
    for col in range(n):
        pivot = next((r for r in range(col, n) if M[r][col] != 0), None)
        if pivot is None:
            return None                                  # no pivot in this column: A is singular
        M[col], M[pivot] = M[pivot], M[col]              # interchange
        p = M[col][col]
        M[col] = [x / p for x in M[col]]                 # scale the pivot to 1
        for r in range(n):
            if r != col and M[r][col] != 0:
                f = M[r][col]
                M[r] = [x - f * y for x, y in zip(M[r], M[col])]   # clear the rest of the column
    return [row[n:] for row in M]                        # the right block is A^-1

def show(M):
    return None if M is None else [[str(x) for x in row] for row in M]

K = [[1, 2, 3], [0, 1, 4], [5, 6, 0]]
print(show(inverse_gauss_jordan(K)))     # [['-24', '18', '5'], ['20', '-15', '-4'], ['-5', '4', '1']]
print(show(inverse_gauss_jordan([[1, 2], [2, 4]])))     # None

A = [[2, 1, -1], [-3, -1, 2], [-2, 1, 2]]                # the running example of Chapter 3
print(show(inverse_gauss_jordan(A)))     # [['4', '3', '-1'], ['-2', '-2', '1'], ['5', '4', '-1']]

Key takeaways

• Computing $A^{-1}$ means solving $A\mathbf{x}_j = \mathbf{e}_j$ for every $j$, all at once with $[,A \mid I,]$.
• Row reduce $[,A \mid I,]$: if the left block becomes $I$, the right block is $A^{-1}$.
• If a zero row appears in the left block, $A$ has fewer than $n$ pivots and is singular.
• The method works because the row operations multiply up to $A^{-1}$.
• Inversion costs about $2n^3$ operations, about three times the cost of solving one system.

🧪 Practice

1. Use row reduction to invert $\begin{bmatrix} 1 & 0 & 1 \ 0 & 2 & 1 \ 1 & 1 & 1 \end{bmatrix}$, and verify the answer.
2. Apply the algorithm to $\begin{bmatrix} 1 & 2 & 3 \ 4 & 5 & 6 \ 7 & 8 & 9 \end{bmatrix}$ and identify the step at which it fails.
3. Invert the lower triangular matrix $\begin{bmatrix} 1 & 0 & 0 \ a & 1 & 0 \ b & c & 1 \end{bmatrix}$ by row reduction. Is the inverse lower triangular, as Section 4.2 predicts?
4. Interview: You need only the third column of $A^{-1}$ for a $5000 \times 5000$ matrix $A$, for example to see how every unknown responds to a change in the third input. How would you get it without computing the whole inverse, and how much work do you save? Hint: the third column of $A^{-1}$ solves one specific linear system.

Properties of Inverses

Why it matters. Rules for inverses let you simplify and rearrange matrix equations without computing anything, the way you rearrange $3x = 6$ without thinking. They also show which manipulations are illegal, such as distributing an inverse over a sum.

Theorem. Let $A$ and $B$ be invertible $n \times n$ matrices and $c \ne 0$ a scalar. Then:

Property	Statement
Inverse of an inverse	$(A^{-1})^{-1} = A$
Product (order reverses)	$(AB)^{-1} = B^{-1}A^{-1}$
Several factors	$(A_1A_2 \cdots A_k)^{-1} = A_k^{-1} \cdots A_2^{-1}A_1^{-1}$
Transpose	$(A^T)^{-1} = (A^{-1})^T$
Scalar multiple	$(cA)^{-1} = c^{-1}A^{-1}$
Powers	$(A^k)^{-1} = (A^{-1})^k$, written $A^{-k}$

Each proof consists of checking that the proposed inverse works, which is enough by uniqueness. For the product,

$$(AB)(B^{-1}A^{-1}) = A(BB^{-1})A^{-1} = AIA^{-1} = AA^{-1} = I.$$

For the transpose, transpose both sides of $AA^{-1} = I$: $(A^{-1})^TA^T = I^T = I$.

The order reversal is the same "socks and shoes" principle as for the transpose. If $AB$ means "do $B$, then $A$", undoing it means "undo $A$, then undo $B$".

Products and invertibility. If $A$ and $B$ are square and $AB$ is invertible, then both $A$ and $B$ are invertible. A product of square matrices is invertible exactly when every factor is.

What does not hold.

• $(A + B)^{-1} \ne A^{-1} + B^{-1}$ in general. Even for numbers, $\frac{1}{1 + 1} \ne \frac11 + \frac11$. The sum of two invertible matrices may even be singular: $I + (-I) = O$.
• $(AB)^{-1} \ne A^{-1}B^{-1}$ unless $A$ and $B$ commute.

Cancellation is restored. In Section 4.1, $AB = AC$ did not imply $B = C$. If $A$ is invertible, it does: multiply both sides on the left by $A^{-1}$. In the same way, $BA = CA$ implies $B = C$ when $A$ is invertible, multiplying on the right. The side matters.

Example: solving a matrix equation. Solve $AXB = C$ for $X$, with $A$ and $B$ invertible. Multiply on the left by $A^{-1}$ and on the right by $B^{-1}$:

$$A^{-1}(AXB)B^{-1} = A^{-1}CB^{-1} \quad\Longrightarrow\quad X = A^{-1}CB^{-1}.$$

Writing $X = A^{-1}B^{-1}C$ or $X = CA^{-1}B^{-1}$ would be wrong.

import numpy as np

rng = np.random.default_rng(5)
A, B, C = (rng.standard_normal((3, 3)) for _ in range(3))   # random matrices are invertible almost surely
inv = np.linalg.inv

print(np.allclose(inv(A @ B), inv(B) @ inv(A)))     # True: order reverses
print(np.allclose(inv(A @ B), inv(A) @ inv(B)))     # False
print(np.allclose(inv(A.T), inv(A).T))              # True
print(np.allclose(inv(2.5 * A), inv(A) / 2.5))      # True
print(np.allclose(inv(A + B), inv(A) + inv(B)))     # False: inverses do not distribute over sums

# Solve AXB = C for X, multiplying on the correct sides
X = inv(A) @ C @ inv(B)
print(np.allclose(A @ X @ B, C))                    # True
print(np.allclose(A @ (inv(A) @ inv(B) @ C) @ B, C))   # False: wrong order of factors

Key takeaways

• $(A^{-1})^{-1} = A$, $(A^T)^{-1} = (A^{-1})^T$, and $(cA)^{-1} = c^{-1}A^{-1}$.
• The inverse of a product reverses the order: $(AB)^{-1} = B^{-1}A^{-1}$.
• A product of square matrices is invertible if and only if every factor is.
• There is no rule for $(A + B)^{-1}$; the sum of invertible matrices can be singular.
• Multiplying by $A^{-1}$ cancels $A$ only from the side it is on.

🧪 Practice

1. Simplify $(A^{-1}B)^{-1}(C^{-1}A)^{-1}$, assuming all matrices are invertible and of the same size.
2. Solve for $X$, assuming the needed inverses exist: (a) $XA = B$, (b) $AX + B = C$, (c) $A^{-1}XA = D$.
3. Prove that if $A$ is invertible and symmetric, then $A^{-1}$ is symmetric.
4. Interview: Show that if $A$, $B$, and $A + B$ are invertible, then $A^{-1} + B^{-1}$ is invertible, and find its inverse. Hint: factor $A^{-1} + B^{-1} = A^{-1}(B + A)B^{-1}$, then use the product rule.

Solving Systems with Inverses

Why it matters. If $A$ is invertible, the system $A\mathbf{x} = \mathbf{b}$ has the explicit solution $\mathbf{x} = A^{-1}\mathbf{b}$. This formula is extremely useful for reasoning: it proves that a solution exists and is unique, and it shows how the solution depends on $\mathbf{b}$. As a numerical method, however, it is almost always the wrong choice, and knowing why is a mark of practical competence.

Theorem. If $A$ is an invertible $n \times n$ matrix, then for every $\mathbf{b} \in \mathbb{R}^n$ the system $A\mathbf{x} = \mathbf{b}$ has exactly one solution, namely $\mathbf{x} = A^{-1}\mathbf{b}$.

Proof. Existence: $A(A^{-1}\mathbf{b}) = (AA^{-1})\mathbf{b} = \mathbf{b}$. Uniqueness: if $A\mathbf{x} = \mathbf{b}$, multiply on the left by $A^{-1}$ to get $\mathbf{x} = A^{-1}\mathbf{b}$. $\blacksquare$

Example. The running example of Chapter 3, $2x + y - z = 8$, $-3x - y + 2z = -11$, $-2x + y + 2z = -3$, has coefficient matrix $A$ with

$$A^{-1} = \begin{bmatrix} 4 & 3 & -1 \ -2 & -2 & 1 \ 5 & 4 & -1 \end{bmatrix}, \qquad A^{-1}\mathbf{b} = \begin{bmatrix} 32 - 33 + 3 \ -16 + 22 - 3 \ 40 - 44 + 3 \end{bmatrix} = \begin{bmatrix} 2 \ 3 \ -1 \end{bmatrix}.$$

Reading the inverse. Since $\mathbf{x} = A^{-1}\mathbf{b} = b_1\operatorname{col}_1(A^{-1}) + \cdots + b_n\operatorname{col}n(A^{-1})$, column $j$ of $A^{-1}$ is the response of the solution to a unit increase in $b_j$, and the entry $(A^{-1}){ij}$ measures how sensitive $x_i$ is to $b_j$. In the Leontief model of Section 3.5, column $j$ of $(I - C)^{-1}$ lists the total production every sector needs to deliver one extra dollar of sector $j$'s product.

Why not to compute $A^{-1}$ in practice.

Approach	Cost for one $\mathbf{b}$	Accuracy
solve(A, b) (elimination, or LU)	about $\tfrac23 n^3$	Backward stable: tiny residual
inv(A) @ b	about $2n^3 + 2n^2$	Often worse, especially for ill-conditioned $A$
LU once, then many solves (Section 4.4)	$\tfrac23 n^3$ once, then about $2n^2$ per $\mathbf{b}$	Same as solve

Computing the inverse does about three times as much work and introduces extra rounding error. The rule of thumb used by numerical analysts: when you see $A^{-1}\mathbf{b}$ in a formula, read it as "solve $A\mathbf{x} = \mathbf{b}$". The inverse itself is worth computing only when its entries are needed, for example when an inverse covariance (precision) matrix is the quantity of interest, or for small matrices of fixed size such as $3 \times 3$ or $4 \times 4$ transforms in graphics.

import numpy as np
from scipy.linalg import hilbert

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
b = np.array([8., -11, -3])
print(np.linalg.inv(A).round(10))      # [[ 4.  3. -1.] [-2. -2.  1.] [ 5.  4. -1.]]
print(np.linalg.inv(A) @ b)            # approximately [ 2.  3. -1.]
print(np.linalg.solve(A, b))           # approximately [ 2.  3. -1.]: same answer, preferred method

# Column j of A^-1 is the response of x to a unit increase in b_j
e2 = np.array([0., 1, 0])
print(np.linalg.solve(A, b + e2) - np.linalg.solve(A, b))   # approximately [ 3. -2.  4.] = column 2 of A^-1

# On an ill-conditioned matrix, the explicit inverse loses much more accuracy
H = hilbert(10)                        # a famously ill-conditioned 10 x 10 matrix
x_true = np.ones(10)
bh = H @ x_true
r_solve = np.linalg.norm(H @ np.linalg.solve(H, bh) - bh)
r_inv   = np.linalg.norm(H @ (np.linalg.inv(H) @ bh) - bh)
print(f"{r_solve:.0e}  {r_inv:.0e}")   # roughly 6e-16 versus 5e-05 (exact values vary by platform)

Key takeaways

• If $A$ is invertible, $A\mathbf{x} = \mathbf{b}$ has exactly one solution for every $\mathbf{b}$: $\mathbf{x} = A^{-1}\mathbf{b}$.
• Column $j$ of $A^{-1}$ describes how the solution responds to a unit change in $b_j$.
• Numerically, solve the system directly instead of forming $A^{-1}$: it is about three times cheaper and more accurate.
• Compute the inverse explicitly only when you need its entries, or for small fixed-size matrices.

🧪 Practice

1. Use $A^{-1}$ for $A = \begin{bmatrix} 2 & 1 \ 5 & 3 \end{bmatrix}$ to solve $A\mathbf{x} = \mathbf{b}$ for $\mathbf{b} = (1, 0)$, $(0, 1)$, and $(3, 7)$.
2. Using the inverse of the running example, find the solution when the right-hand side changes to $(8, -11, 0)$, without solving the system again.
3. Show that if $A$ is invertible, the only solution of $A\mathbf{x} = \mathbf{0}$ is $\mathbf{x} = \mathbf{0}$. Conclude that a square matrix whose homogeneous system has a nontrivial solution cannot be invertible.
4. Interview: A paper writes the ridge regression estimate as $\hat{\mathbf{w}} = (X^TX + \lambda I)^{-1}X^T\mathbf{y}$, and a junior engineer implements it literally with np.linalg.inv. How would you rewrite the computation, and why? Hint: the formula asks for the solution of one linear system with a symmetric matrix; consider solve, or a Cholesky factorization (Section 10.3).

Singular vs. Nonsingular Matrices

Why it matters. Invertibility can be recognized in many apparently different ways: by pivots, by solutions of systems, by columns, by determinants, by eigenvalues. The Invertible Matrix Theorem says that, for square matrices, all of these conditions are the same condition. It ties together everything in the first four chapters, and each later chapter adds new equivalent statements to it.

The Invertible Matrix Theorem. For an $n \times n$ matrix $A$, the following statements are equivalent. Either all of them are true or all of them are false.

	Statement	Viewpoint
(a)	$A$ is invertible	Algebra
(b)	$A$ is row equivalent to $I_n$	Elimination
(c)	$A$ has $n$ pivot positions, i.e. $\operatorname{rank} A = n$	Elimination
(d)	$A\mathbf{x} = \mathbf{0}$ has only the trivial solution	Homogeneous systems
(e)	$A\mathbf{x} = \mathbf{b}$ has at least one solution for every $\mathbf{b}$	Existence
(f)	$A\mathbf{x} = \mathbf{b}$ has exactly one solution for every $\mathbf{b}$	Uniqueness
(g)	The columns of $A$ are linearly independent	Columns (Section 6.3)
(h)	The columns of $A$ span $\mathbb{R}^n$	Columns
(i)	There is an $n \times n$ $C$ with $CA = I$	One-sided inverse
(j)	There is an $n \times n$ $D$ with $AD = I$	One-sided inverse
(k)	$A^T$ is invertible	Rows
(l)	$A$ is a product of elementary matrices	Section 4.4
(m)	$\det A \ne 0$	Chapter 5
(n)	$0$ is not an eigenvalue of $A$	Chapter 8

Why it holds (the core). The key is that $A$ is square. With $n$ rows and $n$ columns:

• $n$ pivots means a pivot in every column, so there are no free variables and solutions are unique: (c) gives (d).
• $n$ pivots also means a pivot in every row, so no row $[,0 \cdots 0 \mid c,]$ can appear and every system is consistent: (c) gives (e).
• Conversely, fewer than $n$ pivots leaves both a column without a pivot (a free variable, so (d) fails) and a row without a pivot (so (e) fails for some $\mathbf{b}$).

For a square matrix, "every system has a solution" and "solutions are unique" therefore stand or fall together. This is false for non-square matrices: a $2 \times 3$ matrix can have the first property but never the second.

Geometric view. A nonsingular matrix maps $\mathbb{R}^n$ onto all of $\mathbb{R}^n$, reshaping space without collapsing any direction. A singular matrix squashes at least one direction to zero, so the outputs fill only a lower-dimensional subspace (the column space, Section 6.6), and many inputs share the same output.

     nonsingular 2 x 2                        singular 2 x 2

     grid --> slanted grid                    grid --> points on a single line
     . . .         .  .  .                    . . .
     . . .  -->   .  .  .                     . . .   -->   . . . . . . .
     . . .       .  .  .                      . . .
     every point has exactly one preimage     some points have no preimage,
                                              others have infinitely many

Nearly singular matrices. In exact arithmetic, a matrix is either singular or not. In floating point, the important question is how close a matrix is to being singular. A matrix whose columns are almost dependent is invertible in theory but behaves badly in practice: solutions are hypersensitive to small changes in the data. The condition number $\kappa(A)$ quantifies this (Section 11.3). A rough rule: with $\kappa(A) \approx 10^k$, expect to lose about $k$ of the 16 significant digits of double precision.

import numpy as np

def classify(A):
    n = A.shape[0]
    r = np.linalg.matrix_rank(A)                        # counts pivots, with a rounding tolerance
    return f"rank {r}/{n}, {'nonsingular' if r == n else 'singular'}, cond {np.linalg.cond(A):.1e}"

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
S = np.array([[1., 2, 3], [2, 4, 6], [1, 0, 1]])        # row 2 = 2 * row 1
print(classify(A))                                       # rank 3/3, nonsingular, cond 4.5e+01
print(classify(S))                                       # rank 2/3, singular, cond about 1e+16 (or inf)

# For singular S, Sx = 0 has a nontrivial solution and some b are unreachable
x = np.array([1., 1, -1])
print(S @ x)                                             # [0. 0. 0.]
b = np.array([1., 0, 0])                                 # needs x1 + 2x2 + 3x3 = 1 and twice that = 0
print(np.linalg.matrix_rank(np.column_stack([S, b])))   # 3 > rank S = 2: Sx = b is inconsistent

# Rounding can hide singularity: this matrix is singular in exact arithmetic
T = np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]])
print(np.linalg.matrix_rank(T))                          # 2
print(np.abs(np.linalg.inv(T)).max() > 1e14)             # True: no error raised, but entries are ~1e16 garbage

Key takeaways

• For a square matrix, invertibility is equivalent to $n$ pivots, to unique solvability of every system, to $A\mathbf{x} = \mathbf{0}$ having only the trivial solution, and to many other conditions.
• For square matrices, existence of solutions for every $\mathbf{b}$ and uniqueness of solutions are the same property.
• Singular matrices collapse at least one direction, so information is lost and cannot be recovered.
• In floating point, ask how close to singular a matrix is; the condition number measures this.

🧪 Practice

1. Use the Invertible Matrix Theorem to decide, with as little computation as possible, whether each matrix is invertible: $\begin{bmatrix} 1 & 2 & 3 \ 0 & 4 & 5 \ 0 & 0 & 6 \end{bmatrix}$, $\begin{bmatrix} 1 & 2 & 3 \ 2 & 4 & 6 \ 7 & 8 & 9 \end{bmatrix}$, $\begin{bmatrix} 3 & 1 \ 6 & 2 \end{bmatrix}$.
2. Let $A$ be a $4 \times 4$ matrix such that $A\mathbf{x} = \mathbf{e}_1$ has two different solutions. What can you conclude about $A\mathbf{x} = \mathbf{e}_2$? Justify with the theorem.
3. Explain why a $3 \times 3$ matrix with a column of zeros, or with two equal rows, must be singular. Which statements of the theorem are easiest to use in each case?
4. Interview: Prove that for square matrices, $AD = I$ implies $DA = I$. Hint: use $AD = I$ to show that $A\mathbf{x} = \mathbf{b}$ has a solution for every $\mathbf{b}$, conclude from the theorem that $A$ is invertible, and multiply $AD = I$ on the left by $A^{-1}$.

<a id="44-elementary-matrices-and-factorizations"></a>

4.4 Elementary Matrices and Factorizations

This subchapter turns Gaussian elimination into matrix algebra: each row operation becomes multiplication by a simple matrix, and the whole elimination becomes a factorization $A = LU$ that can be reused to solve many systems cheaply.

Elementary Matrices

Why it matters. Row operations were introduced in Chapter 3 as procedures performed on equations. Expressing each of them as a matrix lets you reason about elimination with algebra instead of narration: you can compose operations, invert them, and prove theorems about them using the rules of matrix multiplication. Elementary matrices are the building blocks for the factorizations in the rest of this subchapter.

Definition. An elementary matrix is a matrix obtained from the identity $I_n$ by performing exactly one elementary row operation. There are three types, matching the three row operations of Section 3.3. For $n = 3$:

Type	Row operation on $I$	Example	Inverse
Interchange	$R_1 \leftrightarrow R_2$	$\begin{bmatrix} 0 & 1 & 0 \ 1 & 0 & 0 \ 0 & 0 & 1 \end{bmatrix}$	The same matrix
Scaling	$R_3 \leftarrow 5R_3$	$\begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 5 \end{bmatrix}$	Scale by $\tfrac15$: $\operatorname{diag}(1, 1, \tfrac15)$
Replacement	$R_3 \leftarrow R_3 - 4R_1$	$\begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ -4 & 0 & 1 \end{bmatrix}$	$R_3 \leftarrow R_3 + 4R_1$: $\begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ 4 & 0 & 1 \end{bmatrix}$

In general, the replacement $R_i \leftarrow R_i + cR_j$ gives the identity with the extra entry $c$ in position $(i, j)$. The scaling must use $c \ne 0$, exactly as for row operations.

Every elementary matrix is invertible, and its inverse is elementary of the same type. The inverse matrix performs the reverse operation: swap the same two rows back, scale by $1/c$, or subtract the multiple that was added. This is the matrix form of the statement in Section 3.3 that every row operation can be undone.

Structure worth noticing.

• A replacement matrix with $i > j$ (adding a multiple of a higher row to a lower row) is unit lower triangular. These are the only replacements used in the forward phase of elimination, which is why lower triangular matrices appear in the LU factorization.
• An interchange matrix is a permutation matrix: each row and each column contains a single 1.
• A scaling matrix is diagonal.

import numpy as np

def E_swap(n, i, j):
    """Elementary matrix for R_i <-> R_j (0-based indices)."""
    E = np.eye(n)
    E[[i, j]] = E[[j, i]]
    return E

def E_scale(n, i, c):
    """Elementary matrix for R_i <- c R_i, with c != 0."""
    assert c != 0, "scaling by zero is not an elementary operation"
    E = np.eye(n)
    E[i, i] = c
    return E

def E_add(n, i, j, c):
    """Elementary matrix for R_i <- R_i + c R_j, with i != j: the identity plus c in position (i, j)."""
    assert i != j
    E = np.eye(n)
    E[i, j] = c
    return E

print(E_add(3, 2, 0, -4))                                       # [[1 0 0] [0 1 0] [-4 0 1]]

# Each inverse is elementary of the same type, performing the reverse operation
print(np.allclose(np.linalg.inv(E_swap(3, 0, 1)), E_swap(3, 0, 1)))        # True
print(np.allclose(np.linalg.inv(E_scale(3, 2, 5)), E_scale(3, 2, 1 / 5)))  # True
print(np.allclose(np.linalg.inv(E_add(3, 2, 0, -4)), E_add(3, 2, 0, 4)))   # True

Key takeaways

• An elementary matrix is the identity after a single row operation: an interchange, a nonzero scaling, or a replacement.
• The replacement $R_i \leftarrow R_i + cR_j$ is the identity with $c$ in position $(i, j)$.
• Every elementary matrix is invertible, and its inverse is the elementary matrix of the reverse operation.
• Replacements that add multiples of higher rows to lower rows are unit lower triangular.

🧪 Practice

1. Write the $4 \times 4$ elementary matrices for (a) $R_2 \leftrightarrow R_4$, (b) $R_3 \leftarrow -2R_3$, (c) $R_4 \leftarrow R_4 + 7R_2$, and the inverse of each.
2. Which of these are elementary? $\begin{bmatrix} 1 & 0 \ 0 & 0 \end{bmatrix}$, $\ \begin{bmatrix} 1 & 3 \ 0 & 1 \end{bmatrix}$, $\ \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}$, $\ \begin{bmatrix} 2 & 0 \ 0 & 3 \end{bmatrix}$. For each that is not, explain why.
3. Show that the product of the replacement matrices for $R_2 \leftarrow R_2 + aR_1$ and $R_3 \leftarrow R_3 + bR_1$ (both $3 \times 3$) is the identity with $a$ and $b$ in column 1, in either order. Do these two matrices commute?
4. Interview: Is the product of two elementary matrices always elementary? Is the inverse of an elementary matrix always elementary? Give a proof or a counterexample for each. Hint: try multiplying two replacement matrices that act on different pairs of rows.

Row Operations as Matrix Products

Why it matters. This topic connects the procedure of elimination with the algebra of matrices. Once each row operation is a left multiplication, an entire elimination becomes one matrix equation. That equation proves why the $[,A \mid I,]$ method computes inverses, shows that every invertible matrix is a product of elementary matrices, and leads directly to the LU factorization.

Theorem. If $E$ is the elementary matrix obtained by performing a row operation on $I_m$, then for any $m \times n$ matrix $A$, the product $EA$ is the result of performing the same row operation on $A$.

Reason. By the row view of multiplication, $\operatorname{row}_i(EA) = \operatorname{row}_i(E),A$, a combination of the rows of $A$ with weights from row $i$ of $E$. For the replacement $R_3 \leftarrow R_3 - 4R_1$, row 3 of $E$ is $[,-4\ \ 0\ \ 1,]$, so row 3 of $EA$ is $-4\operatorname{row}_1(A) + \operatorname{row}_3(A)$, and the other rows of $E$ copy the other rows of $A$ unchanged. The interchange and scaling cases work the same way.

The recipe is easy to remember: to find the matrix that performs an operation, perform that operation on $I$.

Elimination as a product. The forward phase on the running example from Section 3.3 used three replacements:

$$E_1: R_2 \leftarrow R_2 + \tfrac32 R_1, \qquad E_2: R_3 \leftarrow R_3 + R_1, \qquad E_3: R_3 \leftarrow R_3 - 4R_2.$$

Each multiplies on the left, so the operation performed last is written leftmost:

$$E_3E_2E_1A = \begin{bmatrix} 2 & 1 & -1 \ 0 & \tfrac12 & \tfrac12 \ 0 & 0 & -1 \end{bmatrix} = U.$$

     A  --E1-->  E1 A  --E2-->  E2 E1 A  --E3-->  E3 E2 E1 A = U
                                                  (operations accumulate on the left)

Consequences.

1. Row equivalence is left multiplication. $B$ is row equivalent to $A$ exactly when $B = PA$ for some invertible $P$, a product of elementary matrices.
2. Why $[,A \mid I,] \sim [,I \mid A^{-1},]$ works. If $E_k \cdots E_1A = I$, then $E_k \cdots E_1 = A^{-1}$. Applying the same operations to $[,A \mid I,]$ multiplies both blocks by this product: $[,E_k\cdots E_1A \mid E_k\cdots E_1I,] = [,I \mid A^{-1},]$.
3. Invertible matrices are products of elementary matrices. From $E_k \cdots E_1A = I$, multiply by the inverses: $A = E_1^{-1}E_2^{-1}\cdots E_k^{-1}$, and each $E_i^{-1}$ is elementary. This is statement (l) of the Invertible Matrix Theorem. Every invertible linear transformation is a sequence of swaps, stretches, and shears.

Column operations. Multiplying on the right by an elementary matrix performs the corresponding operation on the columns instead: $AE$ applies to the columns of $A$ the operation that turns the columns of $I$ into the columns of $E$. Column operations do not preserve the solutions of $A\mathbf{x} = \mathbf{b}$. They correspond to a change of variables instead. They are important for determinants (Chapter 5).

import numpy as np
from fractions import Fraction as Fr

def E_add(n, i, j, c):
    """Elementary matrix for R_i <- R_i + c R_j (0-based), with exact Fraction entries."""
    E = [[Fr(int(r == s)) for s in range(n)] for r in range(n)]
    E[i][j] = Fr(c)
    return np.array(E, dtype=object)

def show(M):
    return [[str(x) for x in row] for row in M]

A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=object)
E1 = E_add(3, 1, 0, Fr(3, 2))    # R2 <- R2 + (3/2) R1
E2 = E_add(3, 2, 0, 1)           # R3 <- R3 + R1
E3 = E_add(3, 2, 1, -4)          # R3 <- R3 - 4 R2

U = E3 @ E2 @ E1 @ A             # the first operation is applied first, so it is written rightmost
print(show(U))                   # [['2', '1', '-1'], ['0', '1/2', '1/2'], ['0', '0', '-1']]
print(show(E3 @ E2 @ E1))        # [['1', '0', '0'], ['3/2', '1', '0'], ['-5', '-4', '1']]

# Right multiplication acts on columns: C2 <- C2 - 2 C1 is performed by AE with E[0, 1] = -2
Ac = np.array([[1, 2], [3, 4]])
Ec = np.array([[1, -2], [0, 1]])
print(Ac @ Ec)                   # [[ 1  0] [ 3 -2]]: column 2 minus twice column 1

The accumulated product $E_3E_2E_1$ is unit lower triangular, but its entries ($-5$, $-4$) are not simply the multipliers used. The next topic shows that the inverse product is much tidier.

Key takeaways

• Left-multiplying $A$ by an elementary matrix performs the corresponding row operation on $A$.
• A sequence of row operations is a product $E_k \cdots E_1A$, with the last operation leftmost.
• $B$ is row equivalent to $A$ exactly when $B = PA$ with $P$ invertible.
• Every invertible matrix is a product of elementary matrices, and $E_k\cdots E_1 = A^{-1}$ when the operations reduce $A$ to $I$.
• Right-multiplying by an elementary matrix performs a column operation instead.

🧪 Practice

1. Find elementary matrices $E_1, E_2$ such that $E_2E_1\begin{bmatrix} 1 & 3 \ 2 & 7 \end{bmatrix}$ is upper triangular with $1$s on the diagonal. Then write $\begin{bmatrix} 1 & 3 \ 2 & 7 \end{bmatrix}$ as a product of elementary matrices.
2. A sequence of row operations turns $A$ into $B$. Show that the same sequence turns $AC$ into $BC$ for any matrix $C$ of compatible size. What is the matrix interpretation?
3. Find the single matrix $P$ such that $PA$ swaps rows 1 and 3 of a $3 \times n$ matrix $A$ and then adds twice the new row 1 to row 2.
4. Interview: Explain why Gauss-Jordan elimination on $[,A \mid I,]$ produces $A^{-1}$, using only matrix algebra. Hint: write the whole elimination as a single matrix $P$ multiplying both blocks, and ask what $PA = I$ says about $P$.

LU Decomposition

Why it matters. Gaussian elimination does two things at once: it decides which operations to perform (which depends only on $A$) and it applies them to the right-hand side $\mathbf{b}$. The LU decomposition records the first part once and for all as a factorization $A = LU$. After that, each new right-hand side costs only two cheap triangular solves. LU is the factorization behind most general linear solvers, and it is also the most efficient way to compute determinants of large matrices.

Definition. An LU decomposition (or LU factorization) of a square matrix $A$ is a factorization

$$A = LU,$$

where $L$ is unit lower triangular and $U$ is upper triangular. $U$ is the echelon form produced by the forward phase of elimination, and $L$ stores the multipliers.

Derivation. Suppose the forward phase reaches $U$ using only replacements of the form $R_i \leftarrow R_i - l_{ij}R_j$ with $i > j$, with no interchanges. In matrix form, $E_k \cdots E_1A = U$. Each $E_i$ is unit lower triangular, so

$$A = (E_1^{-1}E_2^{-1}\cdots E_k^{-1}),U = LU,$$

and $L$, a product of unit lower triangular matrices, is unit lower triangular. Something convenient happens in this product: when it is formed in this order, the multipliers do not interact, and each multiplier $l_{ij}$ lands in position $(i, j)$ of $L$. No arithmetic is required to build $L$.

Example. For the running example, the forward phase used $l_{21} = -\tfrac32$ ($R_2 \leftarrow R_2 - (-\tfrac32)R_1$), $l_{31} = -1$, and $l_{32} = 4$. Hence

$$\underbrace{\begin{bmatrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{bmatrix}}{A} = \underbrace{\begin{bmatrix} 1 & 0 & 0 \ -\tfrac32 & 1 & 0 \ -1 & 4 & 1 \end{bmatrix}}{L}\ \underbrace{\begin{bmatrix} 2 & 1 & -1 \ 0 & \tfrac12 & \tfrac12 \ 0 & 0 & -1 \end{bmatrix}}_{U}.$$

Check row 3: $-1\cdot(2, 1, -1) + 4\cdot(0, \tfrac12, \tfrac12) + 1\cdot(0, 0, -1) = (-2, 1, 2)$. Compare $L$ with the tangled product $E_3E_2E_1$ of the previous topic: $L$ is its inverse, and much cleaner.

   During elimination, each multiplier fills the slot it just zeroed out:

   step 1 (column 1)          step 2 (column 2)          result: L and U stored together

   [  2    1   -1 ]           [  2    1    -1  ]         [  2     1     -1  ]    upper part
   [ -3/2  1/2  1/2]          [ -3/2  1/2   1/2]         [ -3/2   1/2   1/2 ]    and diagonal = U
   [ -1    2    1 ]           [ -1    4    -1  ]         [ -1     4     -1  ]    below diagonal = L
     ^ multipliers                  ^ multiplier                                 (unit diagonal implied)

Libraries store $L$ and $U$ in a single array this way, overwriting $A$, so the factorization needs no extra memory.

When does LU exist? Without row interchanges, elimination fails as soon as a zero pivot appears. For example, $\begin{bmatrix} 0 & 1 \ 1 & 1 \end{bmatrix}$ has no LU decomposition. If it had one, $u_{11} = a_{11} = 0$ and $l_{21}u_{11} = a_{21} = 1$, which is impossible. In general, an invertible $A$ has an LU decomposition exactly when every leading principal submatrix (the top-left $k \times k$ corner, for $k = 1, \dots, n$) is invertible, and then $L$ and $U$ are unique. The next topic removes this restriction with row interchanges.

Variant: LDU. Factoring the pivots out of $U$ gives $A = LDU'$ with $D = \operatorname{diag}(u_{11}, \dots, u_{nn})$ and $U'$ unit upper triangular. For a symmetric matrix, $U' = L^T$, so $A = LDL^T$. This is the starting point of the Cholesky factorization (Section 10.3).

Cost. The factorization costs about $\tfrac23 n^3$ operations, the same as the forward phase of elimination, because it is the forward phase.

import numpy as np

def lu_nopivot(A):
    """Doolittle LU without pivoting: A = L U, L unit lower triangular. Fails on a zero pivot."""
    U = np.array(A, dtype=float)
    n = U.shape[0]
    L = np.eye(n)
    for k in range(n - 1):
        if U[k, k] == 0:
            raise ZeroDivisionError(f"zero pivot in column {k + 1}: row interchanges needed")
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]          # the multiplier goes straight into position (i, k)
            U[i, k:] -= L[i, k] * U[k, k:]       # R_i <- R_i - l_ik R_k
    return L, U

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
L, U = lu_nopivot(A)
print(L)                                         # [[ 1.   0.   0. ] [-1.5  1.   0. ] [-1.   4.   1. ]]
print(U)                                         # [[ 2.   1.  -1. ] [ 0.   0.5  0.5] [ 0.   0.  -1. ]]
print(np.allclose(L @ U, A))                     # True

# LDU: factor the pivots out of U
D = np.diag(np.diag(U))
U1 = np.linalg.inv(D) @ U                        # unit upper triangular
print(np.diag(D), np.allclose(L @ D @ U1, A))    # [ 2.   0.5 -1. ] True

try:
    lu_nopivot(np.array([[0., 1], [1, 1]]))
except ZeroDivisionError as e:
    print(e)                                     # zero pivot in column 1: row interchanges needed

Key takeaways

• $A = LU$ with $L$ unit lower triangular (the multipliers) and $U$ upper triangular (the echelon form).
• The multiplier $l_{ij}$ used to clear entry $(i, j)$ goes directly into position $(i, j)$ of $L$.
• LU without interchanges exists exactly when no zero pivot appears; for invertible $A$, this means all leading principal submatrices are invertible.
• The factorization costs about $\tfrac23 n^3$ operations and can overwrite $A$ in place.
• Factoring out the pivots gives $A = LDU'$, and $A = LDL^T$ when $A$ is symmetric.

🧪 Practice

1. Find the LU decomposition of $\begin{bmatrix} 2 & 3 \ 4 & 7 \end{bmatrix}$ and of $\begin{bmatrix} 1 & 2 & 1 \ 2 & 5 & 4 \ 3 & 8 & 9 \end{bmatrix}$, and verify each by multiplying.
2. Write the LDU decomposition of the symmetric matrix $\begin{bmatrix} 4 & 2 \ 2 & 5 \end{bmatrix}$ and confirm that it has the form $LDL^T$.
3. For which values of $c$ does $\begin{bmatrix} c & 1 \ 1 & 1 \end{bmatrix}$ fail to have an LU decomposition? For which values is it singular? Are the two sets the same?
4. Interview: Why does the multiplier $l_{ij}$ land unchanged in position $(i, j)$ of $L = E_1^{-1}\cdots E_k^{-1}$, while $E_k\cdots E_1 = L^{-1}$ has messy entries? Hint: multiply two $3 \times 3$ replacement matrices in both orders, and track whether a row that has already been used as a source is later modified.

PLU Decomposition with Pivoting

Why it matters. Plain LU fails whenever a zero pivot appears, even for perfectly well-behaved invertible matrices, and it is numerically unstable when a pivot is merely small (Section 3.3). Allowing row interchanges fixes both problems. The resulting factorization $PA = LU$ exists for every square matrix and is what libraries such as LAPACK actually compute when you solve a linear system.

Permutation matrices. A permutation matrix $P$ is the identity with its rows reordered. Then $PA$ is $A$ with its rows reordered in the same way. Permutation matrices have useful properties:

• every row and every column contains exactly one 1;
• a product of permutation matrices is a permutation matrix;
• $P^{-1} = P^T$, because the rows of $P$ are distinct standard basis vectors, so $PP^T = I$.

In software, $P$ is stored as an integer vector listing the order of the rows, never as an $n \times n$ matrix.

Theorem. Every $n \times n$ matrix $A$, singular or not, has a factorization

$$PA = LU,$$

where $P$ is a permutation matrix, $L$ is unit lower triangular, and $U$ is upper triangular. Equivalently, $A = P^TLU$.

Partial pivoting. At each column, the algorithm swaps into the pivot position the remaining entry of largest absolute value. This guarantees that every multiplier satisfies $|l_{ij}| \le 1$, which keeps rounding errors from growing (Chapter 12). When rows are swapped, the multipliers already stored in those rows of $L$ are swapped too, so that $L$ stays consistent with the final row order.

Example. Factor $A = \begin{bmatrix} 0 & 2 & 1 \ 1 & 1 & 1 \ 2 & 2 & 4 \end{bmatrix}$ with partial pivoting:

  column 1: largest |entry| is 2 in row 3        column 2: candidates are 0 (row 2) and 2 (row 3)
            swap R1 <-> R3                                 swap R2 <-> R3, multipliers move with the rows

  [ 2  2  4 ]   l21 = 1/2: R2 <- R2 - 1/2 R1    [ 2  2  4 ]               [ 2  2   4 ]
  [ 1  1  1 ]   l31 = 0                   --->  [ 0  0 -1 ]  (l = 1/2)  ---> [ 0  2   1 ]  (l = 0)
  [ 0  2  1 ]                                   [ 0  2  1 ]  (l = 0)       [ 0  0  -1 ]  (l = 1/2)
                                                                           l32 = 0: nothing more to do

  final row order: original rows 3, 1, 2

$$P = \begin{bmatrix} 0 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \end{bmatrix}, \qquad L = \begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ \tfrac12 & 0 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & 2 & 4 \ 0 & 2 & 1 \ 0 & 0 & -1 \end{bmatrix}.$$

Check: row 3 of $LU$ is $\tfrac12(2, 2, 4) + (0, 0, -1) = (1, 1, 1)$, which is row 2 of $A$, as $P$ requires.

Conventions differ. Some references and libraries write $PA = LU$, others $A = PLU$ (then their $P$ is the transpose of ours). SciPy's scipy.linalg.lu returns $P, L, U$ with $A = PLU$. Always check which convention a library uses.

import numpy as np
from scipy.linalg import lu

def plu(A):
    """LU with partial pivoting: returns perm, L, U with A[perm] = L @ U (i.e. PA = LU)."""
    U = np.array(A, dtype=float)
    n = U.shape[0]
    L = np.eye(n)
    perm = np.arange(n)                                  # P stored as a row order, not a matrix
    for k in range(n - 1):
        p = k + np.argmax(np.abs(U[k:, k]))              # largest candidate at or below row k
        if p != k:
            U[[k, p]] = U[[p, k]]                        # swap rows of the partially reduced matrix
            L[[k, p], :k] = L[[p, k], :k]                # swap the multipliers already stored
            perm[[k, p]] = perm[[p, k]]
        if U[k, k] == 0:
            continue                                     # whole column is zero below: singular, skip
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]                  # |L[i, k]| <= 1 thanks to pivoting
            U[i, k:] -= L[i, k] * U[k, k:]
    return perm, L, U

A = np.array([[0., 2, 1], [1, 1, 1], [2, 2, 4]])
perm, L, U = plu(A)
print(perm)                                              # [2 0 1]: original rows 3, 1, 2
print(L)                                                 # [[1.  0.  0. ] [0.  1.  0. ] [0.5 0.  1. ]]
print(U)                                                 # [[ 2.  2.  4.] [ 0.  2.  1.] [ 0.  0. -1.]]
print(np.allclose(A[perm], L @ U))                       # True

P = np.eye(3)[perm]                                      # build the matrix only to check P^-1 = P^T
print(np.allclose(P @ A, L @ U), np.allclose(P @ P.T, np.eye(3)))   # True True

# SciPy uses the convention A = P L U, so its P is the transpose of ours
Ps, Ls, Us = lu(A)
print(np.allclose(Ps, P.T), np.allclose(A, Ps @ Ls @ Us))           # True True

# The matrix with no plain LU factors fine with one interchange
perm2, L2, U2 = plu(np.array([[0., 1], [1, 1]]))
print(perm2, L2.tolist(), U2.tolist())                   # [1 0] [[1.0, 0.0], [0.0, 1.0]] [[1.0, 1.0], [0.0, 1.0]]

Key takeaways

• Every square matrix has a factorization $PA = LU$ with $P$ a permutation matrix.
• Partial pivoting chooses the largest available pivot in each column, which keeps $|l_{ij}| \le 1$ and controls rounding error.
• When rows are swapped mid-elimination, the multipliers already stored in $L$ are swapped too.
• Permutation matrices satisfy $P^{-1} = P^T$ and are stored as index vectors in practice.
• Libraries differ between $PA = LU$ and $A = PLU$; SciPy returns the latter.

🧪 Practice

1. Find a $PA = LU$ factorization of $\begin{bmatrix} 0 & 1 \ 2 & 3 \end{bmatrix}$ and of $\begin{bmatrix} 1 & 2 & 3 \ 2 & 4 & 7 \ 3 & 7 & 11 \end{bmatrix}$, using partial pivoting.
2. Write all six $3 \times 3$ permutation matrices. Which of them equal their own inverse?
3. Apply your PLU procedure to the singular matrix $\begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}$. Where does singularity show up in the factors?
4. Interview: Partial pivoting bounds the multipliers by 1. Can entries of $U$ still grow large during elimination, and why would that matter? Hint: look up the "growth factor" and the classic worst-case matrix with 1 on the diagonal, $-1$ below it, and 1 in the last column; Chapter 12 discusses why such cases are rare in practice.

Solving Systems via LU

Why it matters. Many applications solve $A\mathbf{x} = \mathbf{b}$ for one fixed $A$ and many right-hand sides: a structure analyzed under many load cases, a circuit driven by many input signals, a time-stepping simulation that solves with the same matrix at every step, or an inverse computed column by column. LU does the expensive work once. This answers the interview question from Section 3.2 about 500 right-hand sides.

Two triangular solves. Given $PA = LU$, multiply $A\mathbf{x} = \mathbf{b}$ by $P$ to get $LU\mathbf{x} = P\mathbf{b}$. Write $\mathbf{y} = U\mathbf{x}$ and solve in two steps:

1. Forward substitution: solve $L\mathbf{y} = P\mathbf{b}$ for $\mathbf{y}$.
2. Back substitution: solve $U\mathbf{x} = \mathbf{y}$ for $\mathbf{x}$.

Step 1 replays the forward phase of elimination on $\mathbf{b}$ alone; step 2 is ordinary back substitution. Each step costs about $n^2$ operations.

Example. For the running example, $A = LU$ with no interchanges and $\mathbf{b} = (8, -11, -3)$:

$$L\mathbf{y} = \mathbf{b}: \qquad y_1 = 8, \quad y_2 = -11 + \tfrac32(8) = 1, \quad y_3 = -3 + 1(8) - 4(1) = 1.$$

$$U\mathbf{x} = \mathbf{y}: \qquad -z = 1 \Rightarrow z = -1, \quad \tfrac12 y + \tfrac12(-1) = 1 \Rightarrow y = 3, \quad 2x + 3 - (-1) = 8 \Rightarrow x = 2.$$

The vector $\mathbf{y} = (8, 1, 1)$ is exactly the last column of the augmented matrix after the forward phase in Section 3.3, as it should be.

Cost comparison for $k$ right-hand sides with an $n \times n$ matrix:

Method	Cost	$n = 1000$, $k = 500$
Gaussian elimination from scratch for each $\mathbf{b}$	$k \cdot \tfrac23 n^3$	about $3.3 \times 10^{11}$
Compute $A^{-1}$, then multiply each $\mathbf{b}$	$2n^3 + k \cdot 2n^2$	about $3 \times 10^{9}$
Factor $PA = LU$ once, then two triangular solves per $\mathbf{b}$	$\tfrac23 n^3 + k \cdot 2n^2$	about $1.7 \times 10^{9}$

LU is cheapest and also more accurate than the explicit inverse (Section 4.3). It also gives the determinant almost for free: $\det A = \pm u_{11}u_{22}\cdots u_{nn}$, with the sign set by the number of row interchanges (Chapter 5).

import numpy as np
from scipy.linalg import lu_factor, lu_solve, solve_triangular

A = np.array([[2., 1, -1], [-3, -1, 2], [-2, 1, 2]])
L = np.array([[1., 0, 0], [-1.5, 1, 0], [-1, 4, 1]])
U = np.array([[2., 1, -1], [0, 0.5, 0.5], [0, 0, -1]])
b = np.array([8., -11, -3])

y = solve_triangular(L, b, lower=True, unit_diagonal=True)   # forward substitution
x = solve_triangular(U, y)                                   # back substitution
print(y, x)                                                  # [8. 1. 1.] [ 2.  3. -1.]

# In practice: factor once with pivoting, then reuse for every right-hand side
lu_piv = lu_factor(A)                     # packed L and U in one array, plus pivot indices
B = np.array([[8., 1, 0],                 # three right-hand sides as columns
              [-11, 0, 1],
              [-3, 0, 0]])
X = lu_solve(lu_piv, B)                   # each column solved with two triangular solves
print(np.allclose(A @ X, B))              # True
print(X[:, 0])                            # approximately [ 2.  3. -1.]

# Determinant from U: product of the pivots, sign flipped once per row interchange
lu_packed, piv = lu_piv
swaps = np.sum(piv != np.arange(len(piv)))
print(round((-1) ** swaps * np.prod(np.diag(lu_packed)), 10), round(np.linalg.det(A), 10))   # -1.0 -1.0

Key takeaways

• With $PA = LU$, solve $A\mathbf{x} = \mathbf{b}$ by forward substitution $L\mathbf{y} = P\mathbf{b}$, then back substitution $U\mathbf{x} = \mathbf{y}$.
• Factoring costs about $\tfrac23 n^3$ once; each additional right-hand side costs only about $2n^2$.
• For many right-hand sides, LU beats both repeated elimination and the explicit inverse in cost and accuracy.
• The intermediate vector $\mathbf{y}$ is the right-hand side after the forward phase of elimination.
• The factorization also yields the determinant as the signed product of the pivots.

🧪 Practice

1. Using $L$ and $U$ for the running example, solve $A\mathbf{x} = \mathbf{b}$ for $\mathbf{b} = (1, 0, 0)$. What column of $A^{-1}$ have you found?
2. Use the $PA = LU$ factorization of $\begin{bmatrix} 0 & 2 & 1 \ 1 & 1 & 1 \ 2 & 2 & 4 \end{bmatrix}$ from the previous topic to solve $A\mathbf{x} = (3, 3, 8)$. Remember to permute $\mathbf{b}$ first.
3. A simulation solves $A\mathbf{x}_{t+1} = \mathbf{x}_t$ for 10,000 time steps with a fixed $500 \times 500$ matrix $A$. Estimate the operation counts with and without reusing an LU factorization.
4. Interview: You must compute $\mathbf{c}^TA^{-1}\mathbf{b}$ for a large $A$ and given vectors $\mathbf{b}, \mathbf{c}$. Describe an efficient and accurate way to do it without forming $A^{-1}$. Hint: the expression is $\mathbf{c} \cdot \mathbf{x}$ where $\mathbf{x}$ solves one linear system.

<a id="45-block-matrices"></a>

4.5 Block Matrices

This subchapter treats a large matrix as a small matrix whose entries are themselves matrices, which exposes structure, simplifies proofs, and underlies the fastest algorithms for large problems.

Partitioned Matrix Notation

Why it matters. Large matrices are rarely unstructured. A model with two interacting subsystems produces a matrix with four natural regions: each subsystem's internal behavior, and the two directions of coupling between them. Drawing lines through the matrix to separate these regions, and naming each region, lets you see and use that structure. Block notation is also how high-performance libraries organize computations, by splitting matrices into tiles that fit in fast memory.

Definition. A partitioned (or block) matrix is a matrix divided by horizontal and vertical lines into smaller rectangular submatrices called blocks. A partition of the rows into groups of sizes $m_1, \dots, m_p$ and of the columns into groups of sizes $n_1, \dots, n_q$ produces a $p \times q$ array of blocks, where block $A_{ij}$ has size $m_i \times n_j$:

$$A = \left[\begin{array}{cc|c} 1 & 2 & 0 \ 3 & 4 & 5 \ \hline 6 & 7 & 8 \end{array}\right] = \begin{bmatrix} A_{11} & A_{12} \ A_{21} & A_{22} \end{bmatrix}, \quad A_{11} = \begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix},\ A_{12} = \begin{bmatrix} 0 \ 5 \end{bmatrix},\ A_{21} = \begin{bmatrix} 6 & 7 \end{bmatrix},\ A_{22} = \begin{bmatrix} 8 \end{bmatrix}.$$

All blocks in the same block row have the same number of rows, and all blocks in the same block column have the same number of columns. A matrix can be partitioned in many ways; the right one depends on the problem.

                n1 = 2    n2 = 1
              <-------->  <-->
        ^   [  1    2   |  0  ]       A11: 2 x 2    A12: 2 x 1
   m1=2 |   [  3    4   |  5  ]
        v   [-----------+-----]
   m2=1     [  6    7   |  8  ]       A21: 1 x 2    A22: 1 x 1

Partitions you have already used.

Partition	Notation	Where it appeared
Into columns	$A = [,\mathbf{c}_1\ \ \mathbf{c}_2\ \ \cdots\ \ \mathbf{c}_n,]$	Column view of $A\mathbf{x}$ (Section 3.2)
Into rows	$A = \begin{bmatrix} \operatorname{row}_1(A) \ \vdots \ \operatorname{row}_m(A) \end{bmatrix}$	Row view of $A\mathbf{x}$
Matrix and extra column	$[,A \mid \mathbf{b},]$	Augmented matrix
Two square blocks side by side	$[,A \mid I,]$	Computing inverses (Section 4.3)

Building matrices from blocks. Blocks can also be assembled into a larger matrix. In homogeneous coordinates (Section 2.1), a rotation $R$ ($2 \times 2$) followed by a translation $\mathbf{t}$ ($2 \times 1$) is the single $3 \times 3$ matrix

$$T = \begin{bmatrix} R & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix},$$

where $\mathbf{0}^T$ is a $1 \times 2$ zero row. This compact description says more than the nine numbers it contains.

import numpy as np

A = np.array([[1, 2, 0],
              [3, 4, 5],
              [6, 7, 8]])

# Extract blocks by slicing: rows 0..1 / 2, columns 0..1 / 2
A11, A12 = A[:2, :2], A[:2, 2:]
A21, A22 = A[2:, :2], A[2:, 2:]
print(A11.shape, A12.shape, A21.shape, A22.shape)   # (2, 2) (2, 1) (1, 2) (1, 1)

# Assemble a matrix from blocks; the block sizes must line up
print(np.array_equal(np.block([[A11, A12], [A21, A22]]), A))   # True

# A 2D rigid motion in homogeneous coordinates: [[R, t], [0, 1]]
theta = np.pi / 2
R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
t = np.array([[3.], [1.]])
T = np.block([[R, t], [np.zeros((1, 2)), np.ones((1, 1))]])
p = np.array([1., 0, 1])                            # the point (1, 0) in homogeneous form
print((T @ p).round(10))                            # [3. 2. 1.]: rotated to (0, 1), then moved by (3, 1)

Key takeaways

• A block matrix is a matrix divided by horizontal and vertical lines into submatrices; block $A_{ij}$ has size $m_i \times n_j$.
• Blocks in the same block row share a row count, and blocks in the same block column share a column count.
• Columns, rows, augmented matrices, and $[,A \mid I,]$ are all block partitions.
• Choosing a partition that matches the structure of a problem (subsystems, geometry, data groups) makes that structure visible.

🧪 Practice

1. Partition $\begin{bmatrix} 1 & 0 & 2 & 3 \ 0 & 1 & 4 & 5 \ 0 & 0 & 6 & 7 \ 0 & 0 & 8 & 9 \end{bmatrix}$ into four $2 \times 2$ blocks. Which blocks are special matrices?
2. Write the $4 \times 4$ homogeneous matrix for the 3D transformation "scale every coordinate by 2, then translate by $(1, -1, 5)$" in the block form $\begin{bmatrix} S & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$.
3. A $7 \times 9$ matrix is partitioned with row groups of sizes $3, 4$ and column groups of sizes $2, 5, 2$. List the sizes of all six blocks.
4. Interview: A recommendation system has 10 million users and 1 million items, and the user-item matrix is split into tiles of $10{,}000 \times 10{,}000$ that are distributed across machines. Why tile in two dimensions rather than giving each machine a set of full rows? Hint: compare how much of the other operand each machine must receive when computing a product, for row strips versus square tiles.

Block Addition and Multiplication

Why it matters. Block matrices can be added and multiplied as if their blocks were numbers, as long as the sizes fit and the order inside each product is kept. This makes it possible to compute with structured matrices symbolically, to derive formulas such as block inverses, and to design algorithms that work tile by tile. Most fast matrix-multiplication routines are block algorithms.

Block addition. If $A$ and $B$ have the same size and are partitioned in the same way, then $A + B$ is computed block by block: $(A + B){ij} = A{ij} + B_{ij}$. Scalar multiples work the same way: $(cA){ij} = cA{ij}$.

Block multiplication. Suppose the column partition of $A$ matches the row partition of $B$. Such partitions are called conformable. Then

$$(AB){ij} = \sum{k} A_{ik}B_{kj},$$

exactly like the entry formula, with each product $A_{ik}B_{kj}$ a matrix product. For $2 \times 2$ block matrices:

$$\begin{bmatrix} A_{11} & A_{12} \ A_{21} & A_{22} \end{bmatrix}\begin{bmatrix} B_{11} & B_{12} \ B_{21} & B_{22} \end{bmatrix} = \begin{bmatrix} A_{11}B_{11} + A_{12}B_{21} & A_{11}B_{12} + A_{12}B_{22} \ A_{21}B_{11} + A_{22}B_{21} & A_{21}B_{12} + A_{22}B_{22} \end{bmatrix}.$$

Two cautions:

• Order matters inside each product. Write $A_{11}B_{11}$, never $B_{11}A_{11}$. The blocks do not commute.
• Partitions must conform. The column split of $A$ (for example $2 + 1$) must equal the row split of $B$. The other two splits (rows of $A$, columns of $B$) are free and determine the block layout of $AB$.

          A  (3 x 5)                          B  (5 x 4)

     columns split  2 | 3               rows split   2
                                                    ---
                                                     3
                   \_________________________________/
                     these two splits must be identical

     The row split of A (1 + 2) and the column split of B (3 + 1) are free:
     they give the block layout of AB (block rows 1 + 2, block columns 3 + 1).

Why it works. Entry $(r, s)$ of $AB$ is $\sum_l a_{rl}b_{ls}$. Grouping the indices $l$ according to the column partition of $A$ splits this sum into partial sums, and each partial sum is an entry of one block product $A_{ik}B_{kj}$.

Special cases you already know. The column view $AB = [,A\mathbf{b}_1\ \cdots\ A\mathbf{b}_p,]$ is block multiplication with $A$ unpartitioned and $B$ split into columns. The outer-product view $AB = \sum_k \operatorname{col}_k(A)\operatorname{row}_k(B)$ is $A$ split into columns times $B$ split into rows.

Example: composing rigid motions. With $T_i = \begin{bmatrix} R_i & \mathbf{t}_i \ \mathbf{0}^T & 1 \end{bmatrix}$, block multiplication gives

$$T_1T_2 = \begin{bmatrix} R_1R_2 + \mathbf{t}_1\mathbf{0}^T & R_1\mathbf{t}_2 + \mathbf{t}_1 \ \mathbf{0}^TR_2 + 1\cdot\mathbf{0}^T & \mathbf{0}^T\mathbf{t}_2 + 1 \end{bmatrix} = \begin{bmatrix} R_1R_2 & R_1\mathbf{t}_2 + \mathbf{t}_1 \ \mathbf{0}^T & 1 \end{bmatrix}.$$

Reading the blocks: the combined rotation is $R_1R_2$, and the combined translation is $\mathbf{t}_2$ rotated by $R_1$, plus $\mathbf{t}_1$. The calculation proves this for rotations of any dimension, which an entry-by-entry computation would obscure.

Blocked algorithms. Splitting $n \times n$ matrices into $b \times b$ tiles and multiplying tile by tile does the same arithmetic as the naive method, but each tile fits in the CPU cache and is reused many times while it is there. This reorganization is why optimized BLAS libraries run an order of magnitude faster than a straightforward triple loop.

import numpy as np

rng = np.random.default_rng(6)
A = rng.integers(-3, 4, size=(3, 5))
B = rng.integers(-3, 4, size=(5, 4))

# Column split of A (2 + 3) must match the row split of B (2 + 3); the other splits are free
A11, A12 = A[:1, :2], A[:1, 2:]          # A rows split 1 + 2
A21, A22 = A[1:, :2], A[1:, 2:]
B11, B12 = B[:2, :3], B[:2, 3:]          # B columns split 3 + 1
B21, B22 = B[2:, :3], B[2:, 3:]

AB_blocks = np.block([[A11 @ B11 + A12 @ B21, A11 @ B12 + A12 @ B22],
                      [A21 @ B11 + A22 @ B21, A21 @ B12 + A22 @ B22]])
print(np.array_equal(AB_blocks, A @ B))  # True

def blocked_matmul(A, B, b=64):
    """Tile-by-tile product: the same arithmetic as A @ B, organized for cache reuse."""
    n, m, p = A.shape[0], A.shape[1], B.shape[1]
    C = np.zeros((n, p))
    for i in range(0, n, b):
        for j in range(0, p, b):
            for k in range(0, m, b):     # C_ij += A_ik B_kj, the block entry formula
                C[i:i+b, j:j+b] += A[i:i+b, k:k+b] @ B[k:k+b, j:j+b]
    return C

X, Y = rng.standard_normal((200, 150)), rng.standard_normal((150, 170))
print(np.allclose(blocked_matmul(X, Y), X @ Y))   # True: tile sizes need not divide the dimensions

Key takeaways

• Block matrices partitioned the same way are added block by block.
• If the column partition of $A$ matches the row partition of $B$, then $(AB){ij} = \sum_k A{ik}B_{kj}$, just like the entry formula.
• Keep the order inside every block product; blocks are matrices and do not commute.
• The column, row, and outer-product views of multiplication are special cases of block multiplication.
• Blocked (tiled) multiplication performs the same arithmetic but uses memory caches far better.

🧪 Practice

1. Compute $\begin{bmatrix} I & X \ O & I \end{bmatrix}\begin{bmatrix} I & Y \ O & I \end{bmatrix}$ for square blocks of equal size, and use the result to find $\begin{bmatrix} I & X \ O & I \end{bmatrix}^{-1}$.
2. Let $A = [,A_1\ \ A_2,]$ with $A_1$ of size $m \times k$ and $A_2$ of size $m \times (n - k)$. Express $AA^T$ in terms of $A_1$ and $A_2$.
3. Compute $T^2$ for $T = \begin{bmatrix} I_2 & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$ by block multiplication, then find $T^k$. Interpret the result geometrically.
4. Interview: Strassen's algorithm multiplies two $2 \times 2$ block matrices with 7 block multiplications instead of 8. Applied recursively, what complexity does this give for $n \times n$ matrices, and why is it rarely used for moderate $n$? Hint: solve $T(n) = 7T(n/2) + O(n^2)$, then think about extra additions, memory traffic, and numerical stability.

Block Diagonal and Block Triangular Forms

Why it matters. When a matrix has zero blocks in the right places, a large problem splits into smaller ones. A block diagonal matrix describes subsystems that do not interact at all, so each can be solved separately, possibly in parallel. A block triangular matrix describes a one-way dependency: one subsystem can be solved first, and its result fed into the other. Recognizing these forms, sometimes after reordering the unknowns, is one of the most effective ways to speed up a computation.

Block diagonal matrices. A square matrix of the form

$$A = \begin{bmatrix} A_1 & O & \cdots & O \ O & A_2 & \cdots & O \ \vdots & & \ddots & \vdots \ O & O & \cdots & A_k \end{bmatrix} = \operatorname{diag}(A_1, A_2, \dots, A_k),$$

with each $A_i$ square, is block diagonal. Such matrices behave like diagonal matrices, block by block:

Operation	Result
Product	$\operatorname{diag}(A_i)\operatorname{diag}(B_i) = \operatorname{diag}(A_iB_i)$
Power	$A^m = \operatorname{diag}(A_1^m, \dots, A_k^m)$
Inverse	$A$ is invertible exactly when every $A_i$ is, and $A^{-1} = \operatorname{diag}(A_1^{-1}, \dots, A_k^{-1})$
Linear system	$A\mathbf{x} = \mathbf{b}$ splits into $A_i\mathbf{x}_i = \mathbf{b}_i$, independent for each $i$

Example. The adjacency matrix of a graph with two separate components is block diagonal when the vertices of each component are numbered consecutively. Walks never cross between components, and this is visible in $A^k = \operatorname{diag}(A_1^k, A_2^k)$.

Block triangular matrices. A matrix of the form

$$M = \begin{bmatrix} A & B \ O & D \end{bmatrix},$$

with $A$ and $D$ square, is block upper triangular. The system $M\mathbf{x} = \mathbf{b}$, with $\mathbf{x}$ and $\mathbf{b}$ split to match, reads

$$A\mathbf{x}_1 + B\mathbf{x}_2 = \mathbf{b}_1, \qquad D\mathbf{x}_2 = \mathbf{b}_2.$$

It is solved by block back substitution: first $D\mathbf{x}_2 = \mathbf{b}_2$, then $A\mathbf{x}_1 = \mathbf{b}_1 - B\mathbf{x}_2$. The second subsystem does not depend on the first, but the first depends on the second.

   block diagonal: no coupling              block upper triangular: one-way coupling

   [ A1   0  ] [x1]   [b1]                  [ A   B ] [x1]   [b1]
   [ 0    A2 ] [x2] = [b2]                  [ 0   D ] [x2] = [b2]

   solve A1 x1 = b1  and  A2 x2 = b2        solve D x2 = b2 first,
   independently (in parallel)              then A x1 = b1 - B x2

Invertibility and inverse. $M$ is invertible exactly when both $A$ and $D$ are, and then

$$\begin{bmatrix} A & B \ O & D \end{bmatrix}^{-1} = \begin{bmatrix} A^{-1} & -A^{-1}BD^{-1} \ O & D^{-1} \end{bmatrix}.$$

To verify, multiply out with block multiplication: the top-right block is $A(-A^{-1}BD^{-1}) + BD^{-1} = O$. This generalizes the fact that a triangular matrix is invertible exactly when its diagonal entries are nonzero (Section 4.2), with the diagonal blocks taking the role of the diagonal entries. Block lower triangular matrices behave in the same way.

Finding hidden structure. Reordering the unknowns and equations changes $A$ to $PAP^T$ for a permutation $P$, which may reveal block structure. For example, in a sparse system where some unknowns influence others but not vice versa, a suitable ordering produces a block triangular form, and solvers then handle the blocks one at a time.

import numpy as np
from scipy.linalg import block_diag

A1 = np.array([[2., 1], [1, 3]])
A2 = np.array([[4.]])
A = block_diag(A1, A2)                       # [[2 1 0] [1 3 0] [0 0 4]]
print(np.allclose(np.linalg.inv(A), block_diag(np.linalg.inv(A1), np.linalg.inv(A2))))   # True
print(np.allclose(np.linalg.matrix_power(A, 5),
                  block_diag(np.linalg.matrix_power(A1, 5), A2 ** 5)))                   # True

# Block upper triangular system: solve the bottom block first, then the top
Ab = np.array([[3., 1], [2, 2]])
B  = np.array([[1., 0, 2], [0, 1, 1]])
D  = np.array([[1., 2, 0], [0, 1, 0], [1, 0, 1]])
M = np.block([[Ab, B], [np.zeros((3, 2)), D]])
b = np.arange(1., 6.)
x2 = np.linalg.solve(D, b[2:])               # D x2 = b2
x1 = np.linalg.solve(Ab, b[:2] - B @ x2)     # A x1 = b1 - B x2
print(np.allclose(np.concatenate([x1, x2]), np.linalg.solve(M, b)))   # True

# Block inverse formula for block triangular matrices
inv = np.linalg.inv
M_inv = np.block([[inv(Ab), -inv(Ab) @ B @ inv(D)], [np.zeros((3, 2)), inv(D)]])
print(np.allclose(M_inv, inv(M)))            # True

Key takeaways

• Block diagonal matrices act on each block independently: products, powers, inverses, and systems all split block by block.
• A block upper triangular system is solved by block back substitution, bottom block first.
• $\begin{bmatrix} A & B \ O & D \end{bmatrix}$ is invertible exactly when $A$ and $D$ are, with an explicit block inverse.
• Reordering unknowns ($PAP^T$) can reveal block diagonal or block triangular structure and make a problem much cheaper.

🧪 Practice

1. Find the inverse of $\operatorname{diag}!\left(\begin{bmatrix} 1 & 2 \ 3 & 7 \end{bmatrix}, [5], \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}\right)$.
2. Solve $\begin{bmatrix} 2 & 1 & 1 & 0 \ 1 & 1 & 0 & 1 \ 0 & 0 & 1 & 1 \ 0 & 0 & 0 & 2 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 4 \ 3 \ 3 \ 4 \end{bmatrix}$ by block back substitution with $2 \times 2$ blocks.
3. Derive the formula for $\begin{bmatrix} A & O \ C & D \end{bmatrix}^{-1}$ (block lower triangular), assuming $A$ and $D$ are invertible.
4. Interview: A graph has 1 million vertices split into 1,000 connected components of about 1,000 vertices each. You need to solve a linear system with its Laplacian matrix (Section 13.2). How would you exploit the structure, and roughly how much work does it save compared with a dense solve? Hint: order the vertices by component, then compare $1000 \cdot (1000)^3$ with $(10^6)^3$.

Block Inverses and Schur Complements

Why it matters. Eliminating a whole group of unknowns at once, rather than one unknown at a time, leaves a smaller system for the remaining unknowns. The matrix of that smaller system is the Schur complement. It appears in block inverses and determinants, in domain-decomposition solvers for large simulations, in optimization (reduced KKT systems), and in statistics, where the covariance of one group of Gaussian variables given another is exactly a Schur complement.

Block elimination. Let

$$M = \begin{bmatrix} A & B \ C & D \end{bmatrix},$$

with $A$ square and invertible ($k \times k$), and $D$ square ($l \times l$). Imitate one step of Gaussian elimination, with the block $A$ as the pivot: subtract $CA^{-1}$ times the first block row from the second.

$$\begin{bmatrix} I & O \ -CA^{-1} & I \end{bmatrix}\begin{bmatrix} A & B \ C & D \end{bmatrix} = \begin{bmatrix} A & B \ O & D - CA^{-1}B \end{bmatrix}.$$

The matrix

$$S = D - CA^{-1}B$$

is the Schur complement of $A$ in $M$, sometimes written $M/A$. It is what remains of $D$ after the influence of the first group of unknowns has been eliminated.

Scalar intuition. For $1 \times 1$ blocks, $M = \begin{bmatrix} a & b \ c & d \end{bmatrix}$ and $S = d - \frac{cb}{a} = \frac{ad - bc}{a}$. This is the second pivot of Gaussian elimination, and the determinant appears in the numerator. In general, $\det M = \det A \cdot \det S$ (Chapter 5).

Block LDU factorization. Completing the elimination on the right as well gives

$$M = \begin{bmatrix} I & O \ CA^{-1} & I \end{bmatrix}\begin{bmatrix} A & O \ O & S \end{bmatrix}\begin{bmatrix} I & A^{-1}B \ O & I \end{bmatrix}.$$

The outer factors are always invertible (block unit triangular), so when $A$ is invertible, $M$ is invertible exactly when $S$ is.

Block inverse. Inverting the three factors and multiplying in reverse order gives

$$M^{-1} = \begin{bmatrix} A^{-1} + A^{-1}BS^{-1}CA^{-1} & -A^{-1}BS^{-1} \ -S^{-1}CA^{-1} & S^{-1} \end{bmatrix}.$$

The bottom-right block of $M^{-1}$ is $S^{-1}$. For $1 \times 1$ blocks this gives $\frac{a}{ad - bc}$, in agreement with the $2 \times 2$ formula of Section 4.3.

Solving with the Schur complement. For $M\begin{bmatrix} \mathbf{x}_1 \ \mathbf{x}_2 \end{bmatrix} = \begin{bmatrix} \mathbf{b}_1 \ \mathbf{b}_2 \end{bmatrix}$:

1. Eliminate $\mathbf{x}_1$: solve the reduced system $S\mathbf{x}_2 = \mathbf{b}_2 - CA^{-1}\mathbf{b}_1$.
2. Back substitute: $\mathbf{x}_1 = A^{-1}(\mathbf{b}_1 - B\mathbf{x}_2)$.

This pays off when $A$ is easy to invert (diagonal, block diagonal, or already factored) and $S$ is small. In domain decomposition, $A$ describes the interiors of many subdomains, which decouple into independent blocks, and $S$ describes only the unknowns on the interfaces between them.

   M = [ A  B ]     eliminate the x1 block      [ A   B ]      S = D - C A^-1 B
       [ C  D ]   ------------------------->    [ O   S ]      (a smaller system for x2 alone)

Application: conditioning Gaussians. If a random vector $(\mathbf{x}1, \mathbf{x}2)$ is Gaussian with covariance $\begin{bmatrix} \Sigma{11} & \Sigma{12} \ \Sigma_{21} & \Sigma_{22} \end{bmatrix}$, then the covariance of $\mathbf{x}_1$ given an observed value of $\mathbf{x}_2$ is

$$\Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21},$$

the Schur complement of $\Sigma_{22}$. Observing $\mathbf{x}_2$ subtracts the part of the uncertainty about $\mathbf{x}_1$ that $\mathbf{x}_2$ explains. This formula is at the heart of Kalman filters and Gaussian process regression.

import numpy as np

rng = np.random.default_rng(7)
k, l = 3, 2
M = rng.standard_normal((k + l, k + l)) + 5 * np.eye(k + l)   # diagonal shift keeps the blocks invertible
A, B = M[:k, :k], M[:k, k:]
C, D = M[k:, :k], M[k:, k:]
inv = np.linalg.inv

S = D - C @ inv(A) @ B                                        # Schur complement of A in M

# Block elimination produces a zero block below A
E = np.block([[np.eye(k), np.zeros((k, l))], [-C @ inv(A), np.eye(l)]])
print(np.allclose((E @ M)[k:, :k], 0), np.allclose((E @ M)[k:, k:], S))   # True True

# Block inverse formula
Ai, Si = inv(A), inv(S)
M_inv = np.block([[Ai + Ai @ B @ Si @ C @ Ai, -Ai @ B @ Si],
                  [-Si @ C @ Ai,              Si          ]])
print(np.allclose(M_inv, inv(M)))                             # True

# Determinant: det M = det A * det S
print(np.isclose(np.linalg.det(M), np.linalg.det(A) * np.linalg.det(S)))   # True

# Solve M x = b through the smaller system S x2 = b2 - C A^-1 b1
b = rng.standard_normal(k + l)
b1, b2 = b[:k], b[k:]
x2 = np.linalg.solve(S, b2 - C @ np.linalg.solve(A, b1))
x1 = np.linalg.solve(A, b1 - B @ x2)
print(np.allclose(np.concatenate([x1, x2]), np.linalg.solve(M, b)))       # True

# Conditioning a 2D Gaussian: variance of x1 drops once x2 is observed
Sigma = np.array([[4.0, 1.8],                                 # var(x1) = 4, cov = 1.8, var(x2) = 1
                  [1.8, 1.0]])
cond_var = Sigma[0, 0] - Sigma[0, 1] * Sigma[1, 0] / Sigma[1, 1]
print(round(cond_var, 2))                                     # 0.76: down from 4.0

Key takeaways

• Eliminating the first block of unknowns leaves the Schur complement $S = D - CA^{-1}B$ as the matrix of the reduced system.
• $M = \begin{bmatrix} I & O \ CA^{-1} & I \end{bmatrix}\begin{bmatrix} A & O \ O & S \end{bmatrix}\begin{bmatrix} I & A^{-1}B \ O & I \end{bmatrix}$, so with $A$ invertible, $M$ is invertible exactly when $S$ is.
• The bottom-right block of $M^{-1}$ is $S^{-1}$, and $\det M = \det A \cdot \det S$.
• For $1 \times 1$ blocks, the Schur complement is the second pivot $(ad - bc)/a$ of Gaussian elimination.
• Schur complements give conditional covariances of Gaussians and power domain-decomposition solvers.

🧪 Practice

1. For $M = \begin{bmatrix} 2 & 0 & 1 \ 0 & 2 & 1 \ 1 & 1 & 3 \end{bmatrix}$ with $A$ the top-left $2 \times 2$ block, compute the Schur complement $S$ and use it to find the bottom-right entry of $M^{-1}$.
2. Assuming $D$ is invertible, derive the Schur complement of $D$ in $M$, namely $A - BD^{-1}C$, by eliminating the second block of unknowns instead of the first.
3. Verify the block LDU factorization of $M$ by multiplying the three factors with block multiplication.
4. Interview: Two sensors measure a quantity, and their readings $(x_1, x_2)$ are jointly Gaussian with variances 9 and 4 and covariance 5. After observing $x_2$, what is the remaining variance of $x_1$? What does the answer become if the covariance is 0, and why does that make sense? Hint: the conditional variance is the Schur complement $\Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$.

<a id="5-determinants"></a>