Numerical Thinking and Matrix Computation

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

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,

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:

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. $AB=\left[\begin{array}{ccc}A{b}_1& \cdots & A{b}_p\end{array}\right]$ — 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. $AB={\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.

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.² 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.

The cost of computation

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

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.

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

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

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.

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.

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.

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})$.⁴

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.

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.

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 $\Vert Qx{\Vert }_2^2=(Qx{)}^{\ast }(Qx)=x^{\ast }{Q}^{\ast }Qx=x^{\ast }x=\Vert x{\Vert }_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})$.

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 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.