Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Kronecker product — where it appears

The matrix formed by replacing every entry of one matrix by that entry times the whole of another, so its row and column indices are pairs. Its singular values are all the products of the two sets, and its condition number is therefore the product of theirs.

Named by 11 essays across 3 fields — each of them below, with the objects they name alongside it.

03672108144180216024681012eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix216numbers that describe it108λ smallest0.59λ largest11worst |computed − exact|7.1·10⁻¹³the matrix is never neededand neither is its decomposition

An index that is a pair

A discretisation on a two-dimensional grid of n points a side has n² unknowns and a matrix with n⁴ entries — 10⁸ at n = 100. What that matrix is instead is two Kronecker products of an n × n matrix, which is 2n² numbers, and nothing has been approximated: assembling it was the mistake.

tensor · Kronecker 024681012141610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹index of the singular value at the cutσ ⁄ σ₁eight digits10^-4: 5 Kronecker terms10^-8: 7 Kronecker terms10^-12: 8 Kronecker termsnot closed, and nearly closedrank at the cut8a Kronecker product's1terms at 10⁻⁴5terms at 10⁻⁸7terms a decade0.5the inverse leaves the formatby half a term a decade

A solve that is d decompositions

A Kronecker sum is closed under nothing useful — its inverse is not a Kronecker sum and no factorisation of it is one. What it has instead is eigenvectors that are Kronecker products, so a solve with 1,728 unknowns takes one decomposition of a 12 × 12 matrix and nothing else.

tensor · Kronecker 01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not

An equation whose unknown is a matrix

AX + XB = C is linear in X, so it has a coefficient matrix, and writing it down is the obvious thing to do. At n = 100 that matrix has a hundred million entries for a problem with ten thousand unknowns, and the algorithm everybody uses instead never forms it. Its conditioning is not the eigenvalue gap either, which is the number a reader is invited to consult.

structure · Matrix equation 110¹10²10⁴10⁶10⁸10¹⁰10¹²n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis6 decompositions of an n × nunknowns1.6·10⁴multiplications9.4·10⁵dense factorisation2.5·10¹²fitted exponent7‖Ax − b‖ ⁄ ‖b‖1.8·10⁻¹⁵nothing of size n^dis ever factorised

Five indices are cheaper than two

The same 4,096 unknowns cost 1.049·10⁶ multiplications indexed as a 64 × 64 grid and 1.966·10⁵ indexed as six axes of four. The dense factorisation that ignores the indexing costs 4.581·10¹⁰ at every one of them, and the residual improves in the same direction as the cost.

tensor · Kronecker 357911110¹grid side meigenvalue of C⁻¹Awithin ½ of onethe cluster grows like the sidem = 4: inside of 169m = 6: inside of 3611m = 8: inside of 6413m = 10: inside of 10017the cluster grows with the sideand the spectrum with the area

Four orders of conditioning, and four steps

On a 10×10 grid the two-dimensional kernel's condition number runs from 62 at ρ = 0.5 to 818,561 at ρ = 0.98. The preconditioned step count over the same range runs 18, 21, 21, 22, 21, 19, 18, and the count of eigenvalues the preconditioner actually brings within half a unit of one does not move at all — it is 9, 11, 13, 17 at every correlation the figure will draw.

structure · Toeplitz A6×6B6×6I ⊗ A + Bᵀ ⊗ I36×36one equation, two objectsentries in A and B72entries in the coefficient matrix1296two routes, relative gap1.7·10⁻¹⁶‖AX + XB − C‖/‖C‖1.4·10⁻¹⁶the small squares are the problemand the large one is the notation

The elimination the matrix does not need

The Kronecker form of AX + XB = C is dismissed with a hundred million entries and (2/3)n⁶ operations. Both price an elimination, and after the reduction both routes take, the matrix has exactly n³ nonzeros, none of them above the block diagonal, and nothing left to eliminate.

structure · Matrix equation change in the preconditioned countseparable, circulant, 0.5 → 0.980isotropic, circulant, 0.5 → 0.98190.50.60.70.80.910306090120150180correlation ρconjugate gradient stepsSeparable, no preconditionerSeparable, circulantIsotropic, no preconditionerIsotropic, circulantdashed: no preconditionerthe flat count belonged to the separable kernel

The staircase a separable kernel builds

A block-circulant preconditioner took 18, 21, 22, 21, 19 and 18 steps on a 10 × 10 Toeplitz-block-Toeplitz system as the correlation rose from 0.5 to 0.98 and the unpreconditioned count rose from 43 to 178 — the parameter that makes the problem hard was the one the solver did not notice. That kernel factorises. The isotropic kernel with the same correlation along each axis does not, and on it the preconditioned count climbs 17, 22, 27, 29, 33, 36, while the preconditioner buys a factor of 1.4 where it bought ten. The reason is the spectrum's shape: a product of two one-dimensional spectra is a staircase of ten treads, and a kernel that is not a product gives a ramp of thirty-eight.

structure · Toeplitz steps, ρ = 0.9separable, Chan, side 1628separable, symbol, side 168isotropic, Chan, side 1641isotropic, symbol, side 163546810121416010203040grid sidepreconditioned stepsseparable, Chanseparable, symbolisotropic, Chanisotropic, symbolthe symbol flattens the separable count onlythe isotropic count still grows with the side

The symbol that builds the staircase

A block-circulant preconditioner left the isotropic kernel's spectrum a ramp of thirty-eight clusters where the separable kernel's was a staircase of ten, and the proposal was a circulant sampled from the isotropic kernel's own symbol, to gather the ramp into treads of its own. It gathers the separable kernel instead: six clusters at every grid from 4 to 10 and every correlation, sixty-four of a hundred eigenvalues at exactly one, and eight conjugate-gradient steps at every side from 6 to 16 while Chan's circulant climbs from 18 to 28. On the isotropic kernel it shortens the ramp from 38 clusters to 30 and the step count by about a sixth, and both still grow with the grid. The treads come from a Kronecker product, not from a small difference: on both kernels the matrix and its symbol circulant differ in dozens of directions.

structure · Toeplitz steps, against eight without it6 × 6, a ten-thousandth1510 × 10, a ten-thousandth1316 × 16, a ten-thousandth15010203040fraction of the isotropic kernelconjugate-gradient steps01e-61e-41e-216 × 610 × 1016 × 16left end: the separable kernel alonea millionth is already not nothing

Six steps were six eigenvalues

A circulant sampled from a separable kernel's symbol held conjugate gradients at eight steps on every grid, because the preconditioned spectrum was six values. Add a fraction ε of the isotropic kernel and the prediction was that the six treads would widen, the count would stay near eight while they stayed under one per cent, and then climb towards the isotropic count. A ten-thousandth of isotropy leaves six clusters at the one-per-cent rule and takes the count from 8 to 13–15; a millionth takes it to 11. The eight steps were finite termination on six distinct eigenvalues, not convergence on six clusters, and once the eigenvalues are distinct what a tread costs is set by its width against the solver's tolerance: at a tolerance of 10⁻⁴ a ten-thousandth of isotropy costs two steps. The widths grow in proportion to ε, and summed tread by tread they account for the climb to within four steps. Past a tenth the mixture is harder than either kernel alone — 37 steps on a 16 × 16 grid against 35 for the isotropic kernel — while its condition number is lower than both.

structure · Toeplitz backward error ÷ residualC = 0.5 X, κ(X) 10⁸5.5C = 10⁻² X, κ(X) 10⁸185C = 10⁻⁴ X, κ(X) 10⁸9306C = 10⁻⁶ X, κ(X) 10⁸2.1·10⁵110¹10²10³10⁴10⁵10⁶10⁷κ(X), on a logarithmic scalebackward error ÷ residual10⁰10²10⁴10⁶10⁸C = 0.5 XC = 10⁻² XC = 10⁻⁴ XC = 10⁻⁶ Xdashed: 2√2 μ, the boundthe residual is not the backward error

A backward error the answer does not feel

Both ways of solving AX + XB = C leave a residual at the rounding level, and for a linear system that would be backward stability. For the equation it is not: the smallest perturbation of A, B and C that makes the computed X exact can be far larger. On a family built to make it so — X increasingly ill conditioned, C increasingly small beside AX — the structured backward error of Bartels and Stewart's solution climbs to 2.1·10⁵ times its residual, and the Kronecker system's, whose residual is four times smaller, to 2.0·10⁶, within the bound Higham's μ sets. And the error of the answer does not move. Against the exact X it stays where the size of C puts it, 6·10⁻¹⁵ to 1·10⁻⁹, across eight decades of κ(X). The perturbation that makes the backward error large has to sit on C, along directions X is insensitive to, and so the backward error rises and the answer never feels it.

structure · Matrix equation 468101214160246810nlog₂ of the growth factorW ⊗ I₁W ⊗ I₂W ⊗ I₃look 1 step aheadlook 2 steps aheadlook 3 steps aheadlook 4 steps aheaddashed: 2^(n/d − 1); each family defeats every depth up to its owna fixed depth stretches the exponential

A doubling that arrives two steps late

A pivot rule that looks one elimination ahead chooses Wilkinson's growth row for row, because the doubling it sets up arrives one step later; looking two ahead escapes it. The argument said a matrix whose damage arrives two steps late would defeat the two-step rule the same way. It exists, and it is two Wilkinson matrices interleaved: W ⊗ I₂. Every look-ahead of depth one or two reproduces the greedy rule's growth on it exactly, 4 to 64 from n = 6 to 14, and depth three gives 2. Three copies defeat depth three and yield to four. The growth each depth is held to is 2^(n/d − 1): a fixed depth stretches the exponential and never removes it, at a cost that rises by a factor of about n for every step it looks.

elimination · Elimination

Named alongside it

The objects these essays reach for when they reach for this one.

Condition numberSeparabilityPreconditioningCirculant preconditionerClustered spectrumConjugate gradientsFlop countKronecker sumToeplitz matrixDiscrete laplacianExact ground truthMatrix equation

All concepts