Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Identical algebra, different arithmetic — page 2

Essays 25 to 35 of the 35 on this theme, in the same order.
-47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap The arithmetic underneath

A norm that overflows before it is a norm

The vector of sixteen thousands has a Euclidean norm of 4,000, which fp16 represents exactly. Written as the square root of the sum of squares it returns infinity, because squaring doubles the exponent — and the expression costs half the format's range on the one computation every iterative method performs at every step.

 10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable Elimination, and the swap

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error Elimination, and the swap

The inverse that is never formed

x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.

 10²10³10⁴10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴terms summedrelative error+∞, −∞ 1.01zero 1.00stochastic 0.50nearest 0.47√n against n, fittednearest, fitted exponent0.47stochastic, fitted exponent0.5toward +∞, fitted exponent1twelve seeds averaged at each sizethe slope is the bias, not the precision The arithmetic underneath

The direction the error leans

The size of one rounding error is set by the precision. How ten thousand of them combine is set by something else entirely — the rounding mode — and the fitted exponents are 0.47 for round-to-nearest and 1.01 for round-toward-infinity, on identical data at identical precision.

 10⁻⁴10⁻³10⁻²10⁻¹100.250.50.751anisotropy εsmoothing factor μpoint + fully-line + fullpoint + semi-y⅓, the one-dimensional answertwo routes, three curvesgap between the repairs2.2·10⁻¹⁶scan against closed form2.2·10⁻¹⁶the dashed curve lies on the solid one beneath itone repair, written two ways Iterating, instead of factorising

Smoothing a whole line at once

Solve every grid line in the strong direction exactly rather than sweeping over it, and the smoothing factor goes from 0.9993 back to 0.3340 — which is the one-dimensional answer, on a problem that is not one-dimensional. The repair replaces one ε in the closed form by a one.

 1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit Eigenvalues, singular values, rank

Restarting is a filter

A restart throws away the Ritz values it does not want and begins again from a new starting vector. Written in the eigenbasis, that vector's components have been multiplied by a polynomial with its roots at the discarded values — measured component by component, and agreeing with the polynomial to rounding.

 point · strong in y0.951point · strong in x0.956y-line · strong in y0.037y-line · strong in x0.967semi-y · strong in y0.107semi-y · strong in x0.9661.00 — no convergenceone problem, seen from two sidesy-line, aimed0.037y-line, turned sideways0.97the ratio2631×31 grid, twelve V-cyclesthe direction is in the code, not in the problem Iterating, instead of factorising

Coarsening in one direction only

Leave the smoother alone and halve only the strong direction, and the smoothing factor is 0.3340 — identical to line relaxation's, at every anisotropy and every weight, to twelve digits. The convergence factors are then a factor of three apart, and at 45° both repairs fail outright.

 the grid the operator came froman edge is a coupling the matrix calls strongkeptinterpolatedwhat the entries decidedcoupling ratio, x against y0.001strong couplings across x0strong couplings along y110rows kept or dropped whole1111×11 grid, θ = 0.25the coarse grid, from the matrix alone Iterating, instead of factorising

The coarse grid the matrix chooses

Given a tridiagonal matrix and no information about a grid, the coarsening keeps every other point and derives the weights ½, 1, ½ — the operators the geometric method was handed. Given the anisotropic operator, it discovers semi-coarsening, in the right direction, without a coordinate.

 171319253137434910⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁶|r_nn|8.7·10⁻⁴σ_min3.7·10⁻¹²column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only Eigenvalues, singular values, rank

The cheap rank and what it cannot see

Almost nobody computes singular values to decide a rank. The standard substitute is QR with column pivoting, read off the diagonal of R — and there is a triangular matrix on which the greedy rule makes no interchange at all, has no better column available at any step, and reports a matrix eight orders of magnitude further from singular than it is.

 10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰ Iterating, instead of factorising

One sequence and two recurrences

CGLS and LSQR compute the same iterates — the minimiser over a space is unique, so there is nothing to choose between them in the algebra. At κ = 10⁶ they cost 42 steps and 47. At κ = 10¹⁰ they cost 110 and 209, across four seeds, and the quantity that separates them is the orthogonality of a basis neither of them keeps.

 0246810121410⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸squarings, srelative error against the closed formbest: s = 2the approximant is out of rangethe squarings amplifyboth sides costerror at s = 06·10⁻¹²best, at s = 26.9·10⁻¹⁶error at s = 153.2·10⁻¹²‖A‖ / 2^s at the best1.1the squarings are exact in the algebraand each one doubles the rounding Eigenvalues, singular values, rank

The series that has to be squared back

The Taylor series for the matrix exponential is not wrong — every term is computed correctly — and on Moler and Van Loan's two-by-two its largest term is 5.4 million times the answer it sums to. The method that replaces it scales the matrix down and squares the result back, and both halves of that sentence cost: too few squarings and the approximant is out of range, too many and each one doubles the rounding.

All themes