Numerical Linear Algebra
where the algebra stops being the arithmetic
Series

Circulant — the series

7 essays on one idea, from the one that introduces it to the one that assumes the rest.
  1. 024681012141610⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement2.6·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve4.7·10⁻¹⁶imaginary part of a real spectrum1.2·10⁻¹⁶n = 16, and the whole matrix is 16 numberseigenvectors known in advance

    The matrix that is one row

    A circulant of size 16 is sixteen numbers, has no zero entry anywhere, and hands over its entire spectrum in closed form — the discrete Fourier transform of its first column, exactly. An eigensolver spends a sweep of Jacobi rotations over 256 entries arriving at the same answer, and agrees to 1.2·10⁻¹⁵.

    part 1 · structure
  2. 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10³10⁵σ — how far the periodic wrap is from singularrelative forward error10⁻¹10⁻²10⁻³10⁻⁴10⁻⁵10⁻⁶10⁻⁷10⁻⁸10⁻⁹10⁻¹⁰10⁻¹¹10⁻¹²periodic correctionκ(C)·uelimination on Tn = 64, median of five answersperiodic correction, σ = 10⁻⁸1.2·10⁻⁴elimination on T, σ = 10⁻⁸1.1·10⁻¹⁴κ(C) at σ = 10⁻⁸4·10⁸κ(T) at σ = 10⁻⁸1712the same matrix T in every columnonly the wrap it is solved through changes

    The circulant the problem did not contain

    A matrix that differs from a circulant in two corner entries can be solved through the circulant, by a transform and a two-by-two correction, and the cost claim is exact. The accuracy claim is not. On tridiag(−1, 2 + σ, −1), whose condition number stops at 1,712, the correction is wrong by 1.2·10⁻⁴ at σ = 10⁻⁸ while elimination is right to 1.1·10⁻¹⁴ — because the periodic neighbour is singular at σ = 0 and the two-by-two system inherits that. Solved by Cramer's rule, as here, the two amplifications multiply; a later measurement found that a pivoted solve of the same two-by-two system removes the second.

    part 2 · structure
  3. -7-6-5-4-3-210⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²10⁴log₁₀ of the twist's distance from landing, in turnsrelative forward errorsimple zero, landed on at 0simple zero, landed on at 0.3double zero, d = 2n = 64, twist 10⁻⁵ of a turn from landingpair: error, κ(A) = 4.1·10⁶5.4·10⁻¹²single: error, κ(A) = 3.8·10⁶3·10⁻⁸double: error, κ(A) = 4.2·10¹²7515the same condition number of the wrapa different number of modes near zero

    Two near-zeros cost less than one

    Solve a well-conditioned tridiagonal matrix through a nearly singular wrap and the correction's accuracy is not set by how singular the wrap is. At κ = 4·10⁷ one wrap returns the answer to 8.8·10⁻¹¹ — better than κ·u — and another, at κ = 3.8·10⁷, returns it to 3.3·10⁻⁶. The difference is how many of its samples sit near the symbol's zeros. A real wrap lands on a conjugate pair, a rank-two correction absorbs the pair exactly, and its two-by-two system has condition number 1.00 — which mattered because that system was solved by Cramer's rule; solved with pivoting, the single landing costs what the pair costs.

    part 3 · structure
  4. at σ = 10⁻⁸Cramer's rule1.2·10⁻⁴pivoted2.9·10⁻¹⁰elimination on T1.1·10⁻¹⁴cancellation × u1.7·10⁻¹⁰-10-9-8-7-6-5-4-3-210⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1σ, as a power of tenforward errorCramer's rulepivotedcancellation × uelimination on Tleft: σ small, the wrap nearly singularthe pivoted solve follows the cancellation

    The correction lost to its own two-by-two solve

    Solving a band matrix through a circulant and a small correction was measured losing the answer to 10⁻⁴ where elimination kept it to 10⁻¹⁴, and a wrap that landed one sample on a zero was measured costing four orders more than one that landed a pair. Both measurements solved the two-by-two correction system by Cramer's rule. Solved with a row interchange, the corrected solve on the same matrix loses 2.9·10⁻¹⁰ — the cancellation, and nothing multiplied onto it — and the single landing costs what the pair costs. The loss was in the determinant.

    part 4 · structure
  5. n = 64, twist of half a stepnearest sample, Laplacian0.0024nearest sample, squared5.8·10⁻⁶κ(wrap), squared2.8·10⁶10⁻³10⁻²10⁻¹110⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹θsymbolπ/nLaplacian, θ² near 0its square, θ⁴ near 0no twist puts a sample further than π/n from θ = 0and the symbol's order decides what it sees there

    A zero no twist can step around

    Solved through its circulant wrap with the small capacitance system factorised stably, a banded matrix was found to lose only what the cancellation in the correction costs, and a twist of half a step kept that cancellation near one. That holds only while the cancellation is large. With it kept small, the corrected solve's error is a quarter of the wrap's condition number times u on every case measured — and on a symbol with a zero of fourth order, the square of the Laplacian, no twist can keep the wrap well conditioned, because the nearest sample any twist can reach sees (π/n)⁴. At n = 64 the corrected solve is 34 times less accurate than elimination, where the Laplacian's is four.

    part 5 · structure
  6. at n = 128refined ÷ elimination, n = 1281.4unrefined ÷ elimination3310²10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸matrix size nrelative errorκ(wrap)·ucorrectedeliminationafter one stepthe wrap grows as n⁴ and so does the corrected errorone step puts it on elimination's line

    One step past the zero

    The squared Laplacian, solved through its twisted circulant and a rank-four correction, lost 6 to 34 times more accuracy than elimination, because its symbol's fourth-order zero keeps the wrap as ill conditioned as the matrix whatever the twist. One step of iterative refinement — a residual formed with the band itself, a second solve by the same route — takes it to 0.59, 0.68, 0.73 and 1.40 times elimination's error at n = 16, 32, 64 and 128, and the second and third steps go no lower. The loss the zero imposed was the solver's, and a solver's loss is what refinement removes.

    part 6 · structure
  7. size n163264128(L + s)²1 stepκu 1.2·10⁻¹²1 stepκu 1.9·10⁻¹¹1 stepκu 3.1·10⁻¹⁰1 stepκu 4.9·10⁻⁹(L + s)³1 stepκu 1.2·10⁻¹⁰1 stepκu 7.9·10⁻⁹1 stepκu 5.1·10⁻⁷1 stepκu 3.2·10⁻⁵(L + s)⁴1 stepκu 1.3·10⁻⁸1 stepκu 3.3·10⁻⁶2 stepsκu 8.4·10⁻⁴no first solveκu 2.2·10⁻¹(L + s)⁵3 stepsκu 1.3·10⁻⁶3 stepsκu 1.4·10⁻³no first solveκu 1.4·10⁰no first solveκu 1.4·10³(L + s)⁶3 stepsκu 1.3·10⁻⁴no first solveκu 5.6·10⁻¹no first solveκu 2.3·10³no first solveκu 9.5·10⁶κu: the wrap's condition number times the unit roundoffsteps are not a function of κu

    Where one step stops being enough

    One step of iterative refinement took the squared Laplacian's circulant-wrap solve to elimination's accuracy at every size, and the account was that each step multiplies the error by the wrap's condition number times the rounding, so one step suffices while that product is small. Raised to higher powers, the band tests the account and half of it holds: each step does contract by about κ(wrap)·u. The other half fails. The fifth power at sixteen points needs three steps with κ(wrap)·u near 10⁻⁶, where the third power at 128 points needs one with thirty times more, because the first solve starts up to a thousand times further from the answer than κ(wrap)·u says.

    part 7 · structure

All series