---
title: Block Matrices and the LU Factorization
module: Matrix Algebra
moduleNumber: 2
lessonNumber: 3
order: 203
summary: >
  Partitioning a matrix into blocks lets sums, products, and inverses be computed
  block by block, as if the submatrices were scalars. Block structure also
  underlies the LU factorization A = LU, which splits solving Ax = b into two
  fast triangular solves and repays the cost whenever many systems share one
  coefficient matrix.
topics: [Matrix Algebra]
sources:
  - book: Lay
    ref: "Ch. 2 — Matrix Algebra; §2.4 Partitioned Matrices, §2.5 Matrix Factorizations"
---

Treating a matrix as a list of columns already proved useful in defining the
[matrix product](/linear-algebra/matrix-algebra/matrix-operations). Partitioning
generalizes that view: split a matrix with horizontal and vertical rules into a
grid of smaller matrices, or **blocks**, and the algebra of sums, products, and
inverses carries through block by block. Partitioning makes sparse structure
visible and speeds computation, and it produces the LU factorization used by most
linear-system solvers.

## Partitioned matrices

A matrix may be divided by horizontal and vertical lines into submatrices. For
example,

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

where each $A_{ij}$ is a block, for instance
$A_{11} = \begin{bmatrix} 3 & 0 & -1 \\ -5 & 2 & 4 \end{bmatrix}$ and
$A_{23} = \begin{bmatrix} 7 \end{bmatrix}$. When a matrix models a physical system
— a circuit board of interconnected chips, a sparse aerodynamic grid — a natural
partition groups related variables, and the off-diagonal blocks record the
couplings between groups.

**Addition and scaling.** If $A$ and $B$ have the same size and the same
partition, then $A + B$ is formed block by block, each block being the sum of the
corresponding blocks, and $rA$ scales every block by $r$. Both reduce to the
entrywise operations already defined.

## Block multiplication

The row-column rule works on blocks, provided the partitions line up.

> **Definition (Conformable partitions).** Partitions of $A$ and $B$ are
> **conformable** for block multiplication when the column partition of $A$
> matches the row partition of $B$: the blocks of $A$ have as many columns as the
> matching blocks of $B$ have rows.

Given conformable partitions, $AB$ is computed as though the blocks were scalar
entries, each product written with the block from $A$ on the left, since matrix
multiplication does not commute.

$$
% caption: Block multiplication: each block of AB is a row-of-blocks of A times a
% column-of-blocks of B, exactly the row-column rule one level up.
\begin{tikzpicture}[font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  % A (2x2 blocks)
  \draw[black, thick] (0,0) rectangle (2.0,2.0);
  \draw[black] (1.0,0) -- (1.0,2.0);
  \draw[black] (0,1.0) -- (2.0,1.0);
  \fill[acc!14] (0,1.0) rectangle (2.0,2.0);
  \node at (0.5,1.5) {$A_{11}$}; \node at (1.5,1.5) {$A_{12}$};
  \node at (0.5,0.5) {$A_{21}$}; \node at (1.5,0.5) {$A_{22}$};
  \node[black, font=\scriptsize] at (1.0,2.25) {$A$};
  % B (2x1 blocks)
  \draw[black, thick] (3.0,0) rectangle (4.0,2.0);
  \draw[black] (3.0,1.0) -- (4.0,1.0);
  \fill[acc!14] (3.0,0) rectangle (4.0,2.0);
  \node at (3.5,1.5) {$B_1$}; \node at (3.5,0.5) {$B_2$};
  \node[black, font=\scriptsize] at (3.5,2.25) {$B$};
  % result, drawn as a two-block column (box stands in for the brackets)
  \node[anchor=east] at (4.9,1.0) {$AB =$};
  \draw[black, thick] (5.0,0.2) rectangle (8.5,1.8);
  \draw[black] (5.0,1.0) -- (8.5,1.0);
  \node[font=\scriptsize] at (6.75,1.4) {$A_{11}B_1 + A_{12}B_2$};
  \node[font=\scriptsize] at (6.75,0.6) {$A_{21}B_1 + A_{22}B_2$};
\end{tikzpicture}
$$

Block multiplication is the most general way to view a product: the column form of
$Ax$, the column definition of $AB$, the row-column rule, and the rows-of-$AB$
identity are all special partitions. One more partition gives a form used
repeatedly later.

### Column-row expansion

Partition $A$ into its columns and $B$ into its rows. The blocks are conformable
(one column times one row), and the product becomes a sum of
[outer products](/linear-algebra/matrix-algebra/matrix-operations).[^lay-thm10]

> **Theorem (Column-row expansion).** If $A$ is $m \times n$ and $B$ is
> $n \times p$, then
> $$
> AB = \sum_{k=1}^{n} \operatorname{col}_k(A)\,\operatorname{row}_k(B)
> = \operatorname{col}_1(A)\,\operatorname{row}_1(B) + \cdots
> + \operatorname{col}_n(A)\,\operatorname{row}_n(B).
> $$

Each term $\operatorname{col}_k(A)\operatorname{row}_k(B)$ is an $m \times p$ rank-one
matrix, and their sum reconstructs $AB$. The $(i,j)$-entry of the $k$th term is
$a_{ik}b_{kj}$; summing over $k$ recovers the row-column rule. This decomposition
reappears when a large data matrix is stored as two thin factors, and in the
rank-one layers of the
[SVD](/linear-algebra/symmetric-quadratic-svd/singular-value-decomposition).

## Inverses of block matrices

Block structure often makes an inverse easy to write down.

**Block upper triangular.** Suppose
$A = \begin{bmatrix} A_{11} & A_{12} \\ 0 & A_{22} \end{bmatrix}$ is invertible,
with $A_{11}$ of size $p \times p$ and $A_{22}$ of size $q \times q$. Writing
$A^{-1} = B$ in the same partition and setting $AB = I$ gives four block equations
whose solution is

$$
A^{-1} = \begin{bmatrix} A_{11}^{-1} & -A_{11}^{-1}A_{12}A_{22}^{-1} \\ 0 & A_{22}^{-1} \end{bmatrix}.
$$

The derivation also shows $A_{11}$ and $A_{22}$ must themselves be invertible.

**Block diagonal.** A block diagonal matrix — zero blocks off the diagonal of
blocks — is invertible if and only if each diagonal block is invertible, and its
inverse is block diagonal with each block inverted:

$$
\begin{bmatrix} B & 0 \\ 0 & C \end{bmatrix}^{-1}
= \begin{bmatrix} B^{-1} & 0 \\ 0 & C^{-1} \end{bmatrix}.
$$

These identities let a large problem decouple into independent smaller ones — the
computational reason partitioning appears in high-performance libraries and in
algorithms that estimate
[eigenvalues](/linear-algebra/eigenvalues/eigenvectors-and-eigenvalues).

## The LU factorization

The LU factorization splits a matrix into a lower and an upper triangular factor.

> **Definition (LU factorization).** An $m \times n$ matrix $A$ that can be row
> reduced to echelon form without row interchanges factors as $A = LU$, where $U$
> is an $m \times n$ echelon form of $A$ and $L$ is an $m \times m$ **unit lower
> triangular** matrix (lower triangular with $1$s on the diagonal). $L$ is always
> invertible.

$$
% caption: A = LU: a unit lower triangular L (1s on the diagonal) times an upper
% triangular echelon form U.
\begin{tikzpicture}[font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  % A
  \draw[black, thick] (0,0) rectangle (1.6,1.6);
  \node at (0.8,0.8) {$A$};
  \node at (2.1,0.8) {$=$};
  % L (lower triangular)
  \draw[black, thick] (2.6,0) rectangle (4.2,1.6);
  \fill[acc!14] (2.6,0) -- (2.6,1.6) -- (4.2,0) -- cycle;
  \draw[acc] (2.6,1.6) -- (4.2,0);
  \node[black, font=\scriptsize] at (3.05,0.55) {$L$};
  \node[black, font=\scriptsize] at (3.85,1.35) {$0$};
  \node[black, font=\scriptsize] at (3.4,2.0) {unit lower};
  % U (upper triangular)
  \draw[black, thick] (4.7,0) rectangle (6.3,1.6);
  \fill[acc!14] (4.7,1.6) -- (6.3,1.6) -- (6.3,0) -- cycle;
  \draw[acc] (4.7,1.6) -- (6.3,0);
  \node[black, font=\scriptsize] at (5.85,1.05) {$U$};
  \node[black, font=\scriptsize] at (5.05,0.3) {$0$};
  \node[black, font=\scriptsize] at (5.5,2.0) {echelon};
\end{tikzpicture}
$$

### Two triangular solves

With $A = LU$, the equation $Ax = b$ becomes $L(Ux) = b$. Setting $y = Ux$ splits
the solve into two triangular systems:

$$
Ly = b \quad (\text{forward substitution}), \qquad Ux = y \quad (\text{back substitution}).
$$

Each is cheap because a triangular system is solved by sweeping through the
unknowns in order. The value of the factorization shows when many systems
$Ax = b_1, \dots, Ax = b_p$ share the coefficient matrix $A$: factor once, then
each right-hand side costs only two triangular solves. The
[inverse power method](/linear-algebra/eigenvalues/power-method) and repeated
design iterations both fit this pattern.

$$
% caption: Solving Ax = b in two stages once A = LU is known: forward-solve Ly = b
% for y, then back-solve Ux = y for x.
\begin{tikzpicture}[>=stealth, font=\footnotesize,
  st/.style={draw, minimum width=24mm, minimum height=11mm, align=center}]
  \definecolor{acc}{HTML}{4A6FA5}
  \node[st] (b) at (0,0) {$b$};
  \node[st] (y) at (3.6,0) {solve $Ly = b$\\ gives $y$};
  \node[st] (x) at (7.6,0) {solve $Ux = y$\\ gives $x$};
  \draw[->, thick] (b) -- (y) node[midway, above, font=\scriptsize] {forward};
  \draw[->, thick] (y) -- (x) node[midway, above, font=\scriptsize] {back};
\end{tikzpicture}
$$

$$
% caption: The two solves sweep opposite ways: forward substitution runs down the
% lower-triangular L, back substitution runs up the upper-triangular U.
\begin{tikzpicture}[font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[black, thick] (0,0) rectangle (2.0,2.0);
  \draw[black] (0,2.0) -- (2.0,0);
  \node[black] at (0.6,0.55) {$L$};
  \node[black, font=\scriptsize] at (1.0,2.3) {solve $Ly = b$};
  \draw[->, acc, very thick] (2.35,1.9) -- (2.35,0.1);
  \node[black, font=\scriptsize, anchor=west] at (2.45,1.0) {down};
  \begin{scope}[xshift=4.8cm]
    \draw[black, thick] (0,0) rectangle (2.0,2.0);
    \draw[black] (0,2.0) -- (2.0,0);
    \node[black] at (1.4,1.45) {$U$};
    \node[black, font=\scriptsize] at (1.0,2.3) {solve $Ux = y$};
    \draw[->, acc, very thick] (2.35,0.1) -- (2.35,1.9);
    \node[black, font=\scriptsize, anchor=west] at (2.45,1.0) {up};
  \end{scope}
\end{tikzpicture}
$$

> **Worked example.** Suppose $A = LU$ with
> $L = \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 2 & -5 & 1 & 0 \\ -3 & 8 & 3 & 1 \end{bmatrix}$
> and
> $U = \begin{bmatrix} 3 & -7 & -2 & 2 \\ 0 & -2 & -1 & 2 \\ 0 & 0 & -1 & 1 \\ 0 & 0 & 0 & -1 \end{bmatrix}$.
> To solve $Ax = b$ for $b = (-5, 6, -12, 30)$, first run forward substitution on
> $Ly = b$, top to bottom, each row fixing one $y_i$:
>
> $$
> y_1 = -5, \quad
> y_2 = 6 - (-1)y_1 = 1, \quad
> y_3 = -12 - 2y_1 + 5y_2 = 3, \quad
> y_4 = 30 + 3y_1 - 8y_2 - 3y_3 = -2.
> $$
>
> Then back substitution on $Ux = y$, bottom to top:
>
> $$
> x_4 = \frac{-2}{-1} = 2, \quad
> x_3 = \frac{3 - x_4}{-1} = -1, \quad
> x_2 = \frac{1 + x_3 - 2x_4}{-2} = 2, \quad
> x_1 = \frac{-5 + 7x_2 + 2x_3 - 2x_4}{3} = 1.
> $$
>
> So $x = (1, 2, -1, 2)$. Each triangular solve fixes one new unknown per row,
> which is why the pair costs about $2n^2$ flops against the $\tfrac{2}{3}n^3$ of
> the factorization itself.

### Constructing $L$ and $U$

Suppose $A$ reduces to an echelon form $U$ using only row replacements that add a
multiple of one row to a row below it. Those operations are left-multiplications
by unit lower triangular elementary matrices $E_p \cdots E_1 A = U$, so
$A = (E_p \cdots E_1)^{-1} U = LU$ with $L = (E_p \cdots E_1)^{-1}$, itself unit
lower triangular. In practice $L$ is built without forming any $E_k$: it records
the multipliers used during the reduction.

```algorithm
caption: LU factorization by recording row-replacement multipliers
reduce $A$ to an echelon form $U$ using only row replacements
  that add a multiple of one row to a lower row
set every diagonal entry of $L$ to $1$
for each pivot column $j$, from left to right do
  for each row $i$ below the pivot do
    place in $L_{ij}$ the multiple of the pivot row subtracted from row $i$
      (the entry that was cleared, divided by the pivot)
```

> **Worked example.** Factor
> $A = \begin{bmatrix} 2 & 4 & -1 & 5 & -2 \\ -4 & -5 & 3 & -8 & 1 \\ 2 & -5 & -4 & 1 & 8 \\ -6 & 0 & 7 & -3 & 1 \end{bmatrix}$
> using row replacements only. Clear column $1$ below the pivot $2$ with
> multipliers $-2, 1, -3$ (recorded as $L_{21}, L_{31}, L_{41}$):
>
> $$
> A \sim \begin{bmatrix} 2 & 4 & -1 & 5 & -2 \\ 0 & 3 & 1 & 2 & -3 \\ 0 & -9 & -3 & -4 & 10 \\ 0 & 12 & 4 & 12 & -5 \end{bmatrix}.
> $$
>
> Clear column $2$ below the pivot $3$ with multipliers $-3, 4$
> ($L_{32}, L_{42}$), then clear column $4$ below the pivot $2$ with multiplier
> $2$ ($L_{43}$):
>
> $$
> \sim \begin{bmatrix} 2 & 4 & -1 & 5 & -2 \\ 0 & 3 & 1 & 2 & -3 \\ 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 4 & 7 \end{bmatrix}
> \sim \begin{bmatrix} 2 & 4 & -1 & 5 & -2 \\ 0 & 3 & 1 & 2 & -3 \\ 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 5 \end{bmatrix} = U.
> $$
>
> The recorded multipliers fill the lower triangle of $L$:
>
> $$
> L = \begin{bmatrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 1 & -3 & 1 & 0 \\ -3 & 4 & 2 & 1 \end{bmatrix},
> \qquad
> U = \begin{bmatrix} 2 & 4 & -1 & 5 & -2 \\ 0 & 3 & 1 & 2 & -3 \\ 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 5 \end{bmatrix}.
> $$
>
> Clearing the $-4$ below the first pivot used the multiplier $-4/2 = -2 = L_{21}$;
> the others are recorded the same way. A direct multiplication confirms $LU = A$.

### Row interchanges

When a needed pivot is zero — or, in floating point, when **partial pivoting**
selects the largest available pivot for accuracy — row interchanges enter. The
factorization then takes the **permuted** form $PA = LU$, where $P$ is a
[permutation matrix](/linear-algebra/numerical-linear-algebra/lu-and-cholesky)
recording the swaps, and $L$ is permuted lower triangular. The two-stage solve is
unchanged apart from applying $P$ to $b$ first. Most references to "an LU
factorization" allow this permuted case.

## Operation counts

Flop counts explain why LU is preferred to computing $A^{-1}$. For a dense
$n \times n$ matrix with $n$ moderately large:[^lay-flops]

| Task | Flops | Notes |
| --- | --- | --- |
| LU factorization of $A$ | $\approx \tfrac{2}{3}n^3$ | same cost as reducing $[\,A \mid b\,]$ |
| Solve $Ly = b$ then $Ux = y$ | $\approx 2n^2$ | each triangular solve is $\approx n^2$ |
| Compute $A^{-1}$ | $\approx 2n^3$ | three times the LU cost |
| Multiply $A^{-1}b$ | $\approx 2n^2$ | but often less accurate than the LU solve |

Once $A = LU$ is available, each additional right-hand side costs only $2n^2$
flops, against the $2n^3$ needed to build $A^{-1}$ in the first place. When $A$ is
sparse, $L$ and $U$ often stay sparse while $A^{-1}$ fills in, widening the gap
further. The [numerical module](/linear-algebra/numerical-linear-algebra/lu-and-cholesky)
develops pivoting, stability, and the symmetric Cholesky variant $A = R^\top R$.

[^lay-thm10]: **Lay**, _Linear Algebra and Its Applications_, §2.4 — Theorem 10, the column-row expansion $AB = \sum_k \operatorname{col}_k(A)\operatorname{row}_k(B)$, together with block multiplication of conformable partitions and the inverse formulas for block triangular and block diagonal matrices.
[^lay-flops]: **Lay**, §2.5 — the LU factorization, its two-stage triangular solve, the row-replacement construction of $L$, and the numerical note giving $\tfrac{2}{3}n^3$ flops for the factorization against $2n^3$ for the inverse; partial pivoting yields the permuted factorization $PA = LU$.