---
title: Numerical Thinking and Matrix Computation
module: Numerical Linear Algebra
moduleNumber: 8
lessonNumber: 1
order: 801
summary: >
  Numerical analysis builds efficient discrete algorithms for continuous
  problems, and its cost is dominated as much by memory traffic as by
  arithmetic. Block matrix calculus, flop counts, and the BLAS efficiency ratio
  fix the cost model; triangular and unitary matrices are the two computational
  building blocks every factorization rests on.
topics: [Numerical Linear Algebra]
sources:
  - book: Bornemann
    ref: "Part I — Computing with Matrices; §1 What Is Numerical Analysis?, §2 Matrix Calculus"
  - book: Bornemann
    ref: "§4 Execution Times, §5 Triangular Matrices, §6 Unitary Matrices"
draft: false
---

Exact linear algebra assumes perfect arithmetic. Row reduction produces a
solution, the characteristic polynomial has roots, and the Gram–Schmidt process
returns an orthonormal basis, each on the premise that every operation is carried
out without error. A computer honors none of that premise. It stores numbers to
finite precision, spends far more time moving data than adding it, and cannot
evaluate a limit. Numerical linear algebra studies what becomes of the exact
theory once a machine rounds every operation and spends real time on every memory
access.

## Numerical analysis

> **Definition (Numerical analysis).** The construction and analysis of efficient
> discrete algorithms for solving continuous problems that involve large amounts
> of data.[^bor-na] Here _efficient_ means sparing use of the two scarce
> resources, computation time and memory; _continuous_ means the scalars come
> from $\mathbb{R}$ or $\mathbb{C}$, as in analysis.

Numerical analysis exists because the continuous and the discrete do not match.
An analytical result is the value of a limit, and a limit needs infinite time and
infinite memory. To compute anything at all, continuous
quantities must be discretized into something finite. The three tools for that
are **machine numbers** (finite-precision stand-ins for real numbers),
**iteration** (a finite number of improving steps in place of a limit), and
**approximation** (an answer that is deliberately, and controllably, wrong). The
computed result is allowed to differ from the exact one, provided the error stays
bounded and predictable. Numerical linear algebra is where this program is
learned first, because the problems it treats — solving $Ax = b$, factoring a
matrix, finding eigenvalues — sit underneath nearly all of scientific computing.

Throughout, $\mathbb{K}$ denotes either $\mathbb{R}$ or $\mathbb{C}$, and
$\mathbb{K}^{m \times n}$ the space of $m \times n$ matrices over it. A column
vector is an element of $\mathbb{K}^m = \mathbb{K}^{m \times 1}$; a **co-vector**
(row vector) lives in $\mathbb{K}^{1 \times m}$. The **adjoint** $A^\ast$ is the
transpose $A^\top$ when $\mathbb{K} = \mathbb{R}$ and the conjugate transpose when
$\mathbb{K} = \mathbb{C}$; the co-vector $x^\ast$ is the adjoint of the column
vector $x$.

## Thinking in blocks

Manipulating a matrix as an array of _sub-matrices_ rather than scalars is
clearer on paper and maps onto how fast hardware runs, as the timing table below
shows.

Start from the columns. Writing $e_k$ for the $k$-th standard basis vector, the
columns of $A$ are the images of the basis vectors,

$$
a_k = A e_k \quad (k = 1:n),
$$

using the colon shorthand $k = 1:n$ for $k = 1, 2, \dots, n$. Linearity then turns
a matrix–vector product into a **linear combination of columns**, and a
row-covector times a vector into an **inner product**:

$$
A x = \sum_{k=1}^{n} \xi_k\, a_k,
\qquad
y^\ast x = \sum_{k=1}^{n} \eta_k^\ast\, \xi_k,
\qquad
\lVert x \rVert_2 = \sqrt{x^\ast x}.
$$

The matrix product $C = AB$ admits four equivalent readings, and choosing the
right one is what makes an algorithm fast. Writing $A$ by columns $a_k$ / rows
$a_j^\ast$ and $B$ by columns $b_k$:

- **Column-wise.** $A B = \begin{bmatrix} A b_1 & \cdots & A b_p \end{bmatrix}$ —
  each column of the result is $A$ applied to a column of $B$.
- **Row-wise.** The $j$-th row of $AB$ is $a_j^\ast B$.
- **Rank-one sum.** $A B = \sum_{k=1}^{n} a_k\, b_k^\ast$ — a sum of outer
  products, each an $m \times p$ rank-one matrix.
- **Componentwise.** $\gamma_{jl} = \sum_{k=1}^{n} \alpha_{jk}\beta_{kl}$, the
  familiar formula from a first course.

All four compute the same additions and multiplications in a different order. The
componentwise formula, the one usually taught first, is the least useful here.

> **Worked example.** For
>
> $$
> A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix},
> \qquad
> B = \begin{bmatrix} 5 & 6 \\ 7 & 8 \end{bmatrix},
> $$
>
> the rank-one-sum reading is $AB = a_1 b_1^\ast + a_2 b_2^\ast$, with $a_k$ the
> columns of $A$ and $b_k^\ast$ the rows of $B$:
>
> $$
> \begin{bmatrix}1\\3\end{bmatrix}\begin{bmatrix}5 & 6\end{bmatrix}
> + \begin{bmatrix}2\\4\end{bmatrix}\begin{bmatrix}7 & 8\end{bmatrix}
> = \begin{bmatrix}5 & 6\\15 & 18\end{bmatrix}
> + \begin{bmatrix}14 & 16\\28 & 32\end{bmatrix}
> = \begin{bmatrix}19 & 22\\43 & 50\end{bmatrix}.
> $$
>
> The componentwise formula returns the same entries, e.g. $\gamma_{11} = 1\cdot 5 +
> 2\cdot 7 = 19$; the two readings differ only in which products are summed
> first.

> **Lemma (Block multiplication).** Partition $A \in \mathbb{K}^{m \times n}$ into
> a $q \times r$ grid of blocks $A_{jk}$ and $B \in \mathbb{K}^{n \times p}$ into
> an $r \times s$ grid $B_{kl}$, with compatible dimensions. Then $C = AB$ is the
> $q \times s$ block matrix
>
> $$
> C_{jl} = \sum_{k=1}^{r} A_{jk} B_{kl} \quad (j = 1:q,\; l = 1:s).
> $$

One computes with blocks "as if they were scalars," with one caution: block
factors do not commute, so the order $A_{jk} B_{kl}$ matters and cannot be
rearranged the way scalar entries can.[^bor-block] The four product formulas above
are all special cases of this one lemma, recovered by partitioning $A$ and $B$
into single columns, single rows, or single entries.

$$
% caption: Block multiplication treats each $C_{jl}$ as a scalar-like sum of
% products $A_{jk}B_{kl}$ over the shared block index $k$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% A grid (q=2 x r=3)
\draw[black] (0,0) rectangle (2.4,1.6);
\draw[black] (0.8,0)--(0.8,1.6); \draw[black] (1.6,0)--(1.6,1.6);
\draw[black] (0,0.8)--(2.4,0.8);
\node at (0.4,1.2) {$A_{11}$}; \node[acc] at (1.2,1.2) {$A_{12}$}; \node at (2.0,1.2) {$A_{13}$};
\node at (0.4,0.4) {$A_{21}$}; \node at (1.2,0.4) {$A_{22}$}; \node at (2.0,0.4) {$A_{23}$};
\node at (3.0,0.8) {$\times$};
% B grid (r=3 x s=2)
\draw[black] (3.6,-0.4) rectangle (5.2,1.6);
\draw[black] (4.4,-0.4)--(4.4,1.6);
\draw[black] (3.6,0.4)--(5.2,0.4); \draw[black] (3.6,1.0)--(5.2,1.0);
\node at (4.0,1.3) {$B_{11}$}; \node at (4.8,1.3) {$B_{12}$};
\node[acc] at (4.0,0.7) {$B_{21}$}; \node at (4.8,0.7) {$B_{22}$};
\node at (4.0,0.0) {$B_{31}$}; \node at (4.8,0.0) {$B_{32}$};
\node at (5.8,0.8) {$=$};
% C grid (q=2 x s=2)
\draw[black] (6.4,0) rectangle (8.0,1.6);
\draw[black] (7.2,0)--(7.2,1.6); \draw[black] (6.4,0.8)--(8.0,0.8);
\node[acc] at (6.8,1.2) {$C_{11}$}; \node at (7.6,1.2) {$C_{12}$};
\node at (6.8,0.4) {$C_{21}$}; \node at (7.6,0.4) {$C_{22}$};
\node[acc, align=center] at (7.2,-1.1) {$C_{11}=A_{11}B_{11}+A_{12}B_{21}+A_{13}B_{31}$};
\end{tikzpicture}
$$

## The cost of computation

The natural first cost model counts arithmetic. A **flop** is one real
floating-point operation ($+$, $-$, $\cdot$, $/$, $\sqrt{\;}$). Keeping only the
leading order as dimensions grow:

| Operation | Shapes | Flop (leading order) | $m=n=p$ |
| --- | --- | --- | --- |
| Inner product $x^\ast y$ | $x, y \in \mathbb{R}^m$ | $2m$ | $2m$ |
| Outer product $x y^\ast$ | $x \in \mathbb{R}^m,\, y \in \mathbb{R}^n$ | $mn$ | $m^2$ |
| Matrix–vector $Ax$ | $A \in \mathbb{R}^{m\times n}$ | $2mn$ | $2m^2$ |
| Matrix–matrix $AB$ | $A \in \mathbb{R}^{m\times n}, B \in \mathbb{R}^{n\times p}$ | $2mnp$ | $2m^3$ |

If a flop took a fixed time $t_{\text{flop}}$ and nothing else cost anything,
runtime would just be $T_{\text{peak}} = \#\text{flop} \cdot t_{\text{flop}}$, the
**peak performance** ceiling. Real programs miss this ceiling because of data
movement, not arithmetic.

The gap comes from the **memory hierarchy**. A processor has a few very fast, very
small registers and caches, backed by progressively larger and slower memory, down
to disk. Access speed spans two to three orders of magnitude across the levels.

$$
% caption: The memory hierarchy: each level down is larger and roughly an order
% of magnitude slower, so an algorithm's speed depends on how well it reuses data
% already near the top.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[acc, fill=acc!12] (2.6,3.0) -- (3.4,3.0) -- (3.55,2.4) -- (2.45,2.4) -- cycle;
\node[acc] at (3.0,2.7) {registers};
\draw[black, fill=acc!8] (2.45,2.4) -- (3.55,2.4) -- (3.8,1.6) -- (2.2,1.6) -- cycle;
\node at (3.0,2.0) {L1 / L2 / L3 cache};
\draw[black, fill=acc!6] (2.2,1.6) -- (3.8,1.6) -- (4.1,0.7) -- (1.9,0.7) -- cycle;
\node at (3.0,1.15) {RAM};
\draw[black, fill=acc!4] (1.9,0.7) -- (4.1,0.7) -- (4.5,-0.3) -- (1.5,-0.3) -- cycle;
\node at (3.0,0.2) {SSD / disk};
\draw[->, black] (5.2,-0.3) -- (5.2,3.0);
\node[rotate=90, anchor=south] at (5.55,1.35) {faster, smaller};
\draw[->, black] (0.8,3.0) -- (0.8,-0.3);
\node[rotate=90, anchor=north] at (0.45,1.35) {larger, slower};
\end{tikzpicture}
$$

Let $\#\text{iop}$ count the input/output transfers to main memory and
$t_{\text{iop}}$ their cost. A better model of runtime is

$$
T = \#\text{flop}\cdot t_{\text{flop}} + \#\text{iop}\cdot t_{\text{iop}}
   = T_{\text{peak}}\left(1 + \frac{\tau}{q}\right),
\qquad
\tau = \frac{t_{\text{iop}}}{t_{\text{flop}}} \approx 30,
$$

where the machine constant $\tau$ says a memory access costs roughly thirty flops,
and the algorithm-dependent **efficiency ratio**

$$
q = \frac{\#\text{flop}}{\#\text{iop}} = \text{flops per memory transfer}
$$

is the quantity to maximize. An algorithm with $q$ of order one is memory-bound
and idles near the ceiling; one with $q$ growing with the dimension keeps the
arithmetic units busy.

> **Definition (Efficiency ratio).** For an algorithm, $q = \#\text{flop} /
> \#\text{iop}$ is the number of floating-point operations performed per memory
> transfer. Runtime relative to peak is governed by $1 + \tau/q$; large $q$ hides
> the memory-access penalty $\tau$.[^bor-cost]

The three product formulas differ precisely in $q$. An inner or outer product has
$q \approx 1$; a matrix–vector product has $q \approx 2$; a matrix–matrix product
has $q \approx 2m/3$, growing with the dimension. This is why the same arithmetic,
reordered, runs an order of magnitude apart on real hardware.

| Operation | Flop | Memory transfers | $q$ |
| --- | --- | --- | --- |
| Inner product $x^\ast y$ | $2m$ | $2m$ | $1$ |
| Outer product $x y^\ast$ | $m^2$ | $m^2$ | $1$ |
| Matrix–vector $Ax$ | $2m^2$ | $m^2$ | $2$ |
| Matrix–matrix $AB$ | $2m^3$ | $3m^2$ | $2m/3$ |

The **BLAS** (Basic Linear Algebra Subprograms) library packages hardware-tuned
kernels at three levels matched to this table: Level 1 for vector operations (one
loop), Level 2 for matrix–vector operations (two loops), Level 3 for
matrix–matrix operations (three loops). The design rule that follows is to phrase
an algorithm at the highest BLAS level available, so it spends its time in
$O(m^3)$ arithmetic with $O(m^2)$ traffic rather than the reverse.

One more hardware fact leaks into the mathematics. Fortran, and therefore BLAS,
LAPACK, and MATLAB, store a matrix **column by column** (column-major order);
C and Python store it **row by row**. Sweeping a matrix along its stored direction
avoids index arithmetic and cache misses, which is why the column-oriented product
formulas are the default in this subject.

## Triangular matrices

Two families of matrices are the reusable pieces every factorization is built
from. The first is triangular.

> **Definition (Triangular matrix).** $L \in \mathbb{K}^{m \times m}$ is **lower
> triangular** if all entries above the diagonal vanish, and $U$ is **upper
> triangular** if all entries below it vanish. A triangular matrix is **unipotent**
> if every diagonal entry equals $1$.

Two structural facts make them computationally convenient:

- **Determinant.** By Laplace expansion the determinant of a triangular matrix is
  the product of its diagonal entries, $\det = \lambda_1 \cdots \lambda_m$. So a
  triangular matrix is invertible exactly when no diagonal entry is zero.
- **Closure.** Invertible lower (upper) triangular matrices are closed under
  product and inverse; they form a subgroup of the general linear group
  $\mathrm{GL}(m; \mathbb{K})$.[^bor-tri]

The payoff is that a triangular system is solved directly, with no elimination.
For $Lx = b$ with $L$ lower triangular, partition off the last row at each step:
the equation $l_{k-1}^\ast x_{k-1} + \lambda_k \xi_k = \beta_k$ gives $\xi_k$ once
the earlier components are known. Sweeping $k = 1, 2, \dots, m$ is **forward
substitution**; sweeping an upper-triangular system from the bottom up is **back
substitution**.

$$
\xi_k = \frac{\beta_k - l_{k-1}^\ast\, x_{k-1}}{\lambda_k} \quad (k = 1:m).
$$

$$
% caption: Forward substitution fills the solution top to bottom; solving row $k$
% uses only the already-known components above it (shaded) against the lower
% triangle of $L$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% L matrix as lower triangle
\draw[black] (0,0) rectangle (2.4,2.4);
\fill[acc!12] (0,2.4)--(2.4,0)--(0,0)--cycle;
\draw[black] (0,2.4)--(2.4,0);
\node[black] at (0.55,1.85) {$L$};
% highlight solved row k=3 band
\draw[acc, thick] (0,0.8) rectangle (2.4,1.2);
\node[acc] at (3.0,1.0) {row $k$};
% x vector
\draw[black] (4.2,0) rectangle (4.7,2.4);
\fill[acc!25] (4.2,1.2) rectangle (4.7,2.4);
\node[acc] at (4.45,1.8) {known};
\node at (4.45,0.6) {$?$};
\node at (5.1,1.2) {$x$};
% b vector
\node at (5.7,1.2) {$=$};
\draw[black] (6.1,0) rectangle (6.6,2.4);
\node at (6.35,1.2) {$b$};
\draw[->, acc, thick] (0.3,2.1) -- (0.3,0.3);
\node[acc, anchor=west] at (-0.05,-0.35) {sweep top to bottom};
\end{tikzpicture}
$$

```algorithm
caption: $\textsc{ForwardSubstitution}(L, b)$ — solve $Lx = b$ for lower triangular $L$
for $k = 1$ to $m$ do
  $s \gets \beta_k - l_{k-1}^\ast x_{k-1}$ // inner product over the solved part
  $\xi_k \gets s / \lambda_k$
end for
return $x$
```

> **Worked example.** Solve $Lx = b$ by forward substitution with
>
> $$
> L = \begin{bmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ 4 & 1 & 2 \end{bmatrix},
> \qquad
> b = \begin{bmatrix} 4 \\ 11 \\ 19 \end{bmatrix}.
> $$
>
> Row 1 reads $2\xi_1 = 4$, giving $\xi_1 = 2$. Row 2 reads $1\cdot 2 + 3\xi_2 =
> 11$, so $3\xi_2 = 9$ and $\xi_2 = 3$. Row 3 reads $4\cdot 2 + 1\cdot 3 + 2\xi_3 =
> 19$, so $2\xi_3 = 8$ and $\xi_3 = 4$. Each row uses only the components already
> found above it: $x = (2, 3, 4)^\top$.

The inner product at row $k$ costs $2k$ flops, so the whole solve costs
$\sum_{k=1}^m 2k \approx m^2$ flops — the same leading cost, and the same
efficiency ratio $q \approx 2$, as a single triangular matrix–vector product.
Both are standardized Level-2 BLAS routines. Because $\beta_k$ is never needed
after $\xi_k$ is computed, the output can overwrite the input in place, saving a
copy; such **in situ** execution is the BLAS norm.

## Unitary matrices

The second building block preserves geometry rather than exploiting sparsity.

> **Definition (Unitary matrix).** $Q \in \mathbb{K}^{m \times m}$ is **unitary**
> (for $\mathbb{K} = \mathbb{R}$, **orthogonal**) if $Q^{-1} = Q^\ast$,
> equivalently $Q^\ast Q = Q Q^\ast = I$. Its columns then satisfy $q_j^\ast q_l
> = [\,j = l\,]$, i.e. they form an orthonormal basis of $\mathbb{K}^m$.

Three properties matter downstream:

- **Trivial inverse.** Solving $Qx = b$ needs no factorization: $x = Q^\ast b$,
  one matrix–vector product at $2m^2$ flops, twice a triangular solve.
- **Length preservation.** From $\lVert Qx \rVert_2^2 = (Qx)^\ast(Qx) = x^\ast
  Q^\ast Q x = x^\ast x = \lVert x \rVert_2^2$, a unitary map is an isometry: it
  rotates or reflects without stretching. Length preservation keeps errors from
  being amplified, which is why unitary matrices anchor the stable algorithms of
  numerical linear algebra.
- **Group.** Unitary matrices are closed under product and inverse, forming the
  subgroup $U(m)$ ($O(m)$ in the real case) of $\mathrm{GL}(m;\mathbb{K})$.

$$
% caption: A unitary map sends the unit circle to itself and an orthonormal frame
% to another orthonormal frame: lengths and right angles are preserved, only the
% orientation of the frame changes.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% left: standard frame
\draw[black] (-1.3,0) -- (1.3,0);
\draw[black] (0,-1.3) -- (0,1.3);
\draw[black] (0,0) circle (1);
\draw[->, acc, thick] (0,0) -- (1,0); \node[acc, anchor=west] at (1.05,0) {$e_1$};
\draw[->, acc, thick] (0,0) -- (0,1); \node[acc, anchor=south] at (0,1.05) {$e_2$};
% arrow
\draw[->, black, thick] (1.9,0) -- (3.1,0);
\node[anchor=south] at (2.5,0.1) {$Q$};
% right: rotated frame
\begin{scope}[xshift=5cm]
\draw[black] (-1.3,0) -- (1.3,0);
\draw[black] (0,-1.3) -- (0,1.3);
\draw[black] (0,0) circle (1);
\draw[->, acc, thick] (0,0) -- (0.71,0.71); \node[acc, anchor=south west] at (0.71,0.71) {$q_1$};
\draw[->, acc, thick] (0,0) -- (-0.71,0.71); \node[acc, anchor=south east] at (-0.71,0.71) {$q_2$};
\end{scope}
\end{tikzpicture}
$$

A **permutation matrix** $P_\pi$, whose columns are the basis vectors in permuted
order, is the simplest unitary matrix: it merely reorders rows or columns and
satisfies $P_\pi^\ast = P_\pi^{-1}$. Permutations are the bookkeeping of pivoting.

## Building blocks for the factorizations

Triangular and unitary matrices are the two shapes a hard matrix is reduced to. A
triangular factorization ($A = LU$, or Cholesky's $A = R^\ast R$) turns a linear
system into two cheap triangular solves. The [QR
factorization](/linear-algebra/numerical-linear-algebra/qr-and-numerical-least-squares)
reduces a matrix to a unitary times a triangular factor, trading a little extra
cost for the stability that unitary maps guarantee. Each such algorithm is
assembled from block products, forward and back substitution, and orthogonal
transformations, phrased at the highest BLAS level the hardware rewards.

[^bor-na]: **Bornemann**, _Numerical Linear Algebra_, §1 — What Is Numerical Analysis?: the definition in terms of efficient discrete algorithms for continuous problems, and the roles of machine numbers, iteration, and approximation.
[^bor-block]: **Bornemann**, §2.10–§2.17 — Matrix Calculus: the four product formulas (column-wise, row-wise, rank-one sum, componentwise) and the block-multiplication lemma that subsumes them.
[^bor-cost]: **Bornemann**, §4 — Execution Times: the flop table, the runtime model $T = T_{\text{peak}}(1 + \tau/q)$, the efficiency ratio $q$, and the three BLAS levels; §4.8 on column-major versus row-major storage.
[^bor-tri]: **Bornemann**, §5 — Triangular Matrices: the determinant and closure properties, forward and back substitution, and their $m^2$ cost; §6 — Unitary Matrices: $Q^{-1} = Q^\ast$, length preservation, and permutation matrices.