Generator

krylov-exp

One function in the funm library, called 15 times across 13 essays. Below: what it draws at its defaults, what it draws at every value an essay asks for, the 7 claims it put to the test while drawing them, and where it stands against the rule this site is named for.

At its defaults it draws relative error of e^ab by krylov against the number of matrix–vector products, n = 100. The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.

krylov-exp is one function in lib/figures/funm.js — matrix functions — the definition that is not a method, and the vector that was wanted. Everything below came out of it during this build, at arguments taken from the essays rather than invented for this page. A figure here is the figure a reader meets in an essay, and if the generator changes, this page changes with it.

At its defaults

Drawn even though every essay passes arguments — which on this site is every essay, at 100% of placements since the standard pass. A default nothing exercises is a trap for the next essay to call this with none, and this is the page where a default that has drifted from the figures around it becomes visible.

Relative error of e^Ab by Krylov against the number of matrix–vector products, n = 100The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 3·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension100Krylov megaflops0.36dense megaflops2the exponential that is computedis 18×18

The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.

n: 100

The arguments are the ones A parameter that counts steps passes. A value drawn at the generator's defaults instead would be a picture no essay asked for and no assertion has been run against.

Relative error of e^Ab by Krylov against the number of matrix–vector products, n = 100The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 3·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension100Krylov megaflops0.36dense megaflops2the exponential that is computedis 18×18

The exponential of a 100×100 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.026, 8.4·10⁻⁴, 1.5·10⁻⁵, 1.6·10⁻⁷, 1.1·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 3·10⁻¹⁶, at m = 18. That is 0.36 against 2.00 megaflops, on a dense matrix; on a sparse one the ratio is far larger.

n: 160

The arguments are the ones Restarting is a filter passes. A value drawn at the generator's defaults instead would be a picture no essay asked for and no assertion has been run against.

Relative error of e^Ab by Krylov against the number of matrix–vector products, n = 160The exponential of a 160×160 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.027, 8.7·10⁻⁴, 1.5·10⁻⁵, 1.7·10⁻⁷, 1.2·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 2.9·10⁻¹⁶, at m = 18. That is 0.92 against 8.19 megaflops, on a dense matrix; on a sparse one the ratio is far larger.0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in e^Abforming e^A: 2.9·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension160Krylov megaflops0.92dense megaflops8.2the exponential that is computedis 18×18

The exponential of a 160×160 matrix is never formed: m products with A build a Krylov basis, the exponential of the m×m Hessenberg matrix is taken, and one combination of the basis vectors is the answer. The error falls superlinearly — 0.027, 8.7·10⁻⁴, 1.5·10⁻⁵, 1.7·10⁻⁷, 1.2·10⁻⁹ at m = 4, 6, 8, 10, 12 — and crosses the accuracy of the full dense exponential, drawn as the flat line at 2.9·10⁻¹⁶, at m = 18. That is 0.92 against 8.19 megaflops, on a dense matrix; on a sparse one the ratio is far larger.

What it checked while drawing

Every figure above asserted its own claims on the way to being drawn, and a claim that failed would have failed the build rather than drawn a wrong picture. Those assertions used to leave no trace at all: a passing one returned true and the only evidence the figure had checked anything was that nothing crashed. The list below is what they actually said, collected by running this generator with an observer installed — not a description of what it is believed to check.

7 distinct claims across 3 sets of arguments, grouped below by shape — because most of them are one sentence with a different number in it, and how many separate times that sentence was put to the test is the informative part.

a size the dense exponential is still affordable at

a superdiagonal the closed form stays representable at

and the convergence is superlinear rather than a fixed rate

in half the dimension or less

LU is for square matrices

matmul shapes agree

the Krylov route reaches the accuracy of the full exponential inside the sweep

Against the rule

The rule does not apply to it. It factorises nothing, so there is no residual it could be withholding. That is worth stating rather than leaving blank: a site that reported the rule as satisfied by every generator would be counting mostly generators the rule never reached.

Across the library: the rule bites on 92 of 198 generators — 77 print a residual and 15 are exempt with a published reason; 106 factorise nothing. Read from lib/residual-rule.js, which is the same body the gate enforces from, and the gate's last check fails the build if this page and it disagree about any generator.

Where it is called

Changing this generator changes every figure on this list. That is what makes the list worth publishing rather than keeping in a check script.

Methods that were designed apart

A parameter that counts steps

The regularisation field's knob is a positive real number chosen by one of three rules. The iterative field's is an integer nobody called a knob — where to stop. On the same problem the best step is 20 and the best λ is 0.025, and they reach 0.1426 and 0.1406.

Eigenvalues, singular values, rank

An eigenvalue that arrives twice

A matrix with forty distinct eigenvalues, handed to Lanczos for eighty steps, returns twenty-five extra copies of thirteen of them — the largest arriving five times. Every copy is accurate to 1.9·10⁻⁸ relative. No arithmetic error was made, nothing overflowed, and a caller counting eigenvalues gets the wrong multiplicity from a computation in which no individual number is wrong.

Two errors, and whose fault they are

An estimate that can be fooled

Nobody computes a condition number, because forming an inverse costs more than the solve did. Every library estimates it instead, from four or five products with a factorisation already in hand. The estimate is exactly right on four random matrices out of five — and there is a matrix, three distinct entries wide, on which it returns a twentieth of the truth.

Iterating, instead of factorising

An operator with no entries

At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.

Iterating, instead of factorising

An orthogonalisation nobody calls one

Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.

Randomised, and the guarantee that changes kind

Counting what cannot be looked at

The trace is n additions and one of the most expensive quantities in the subject to estimate, because the matrices whose trace is wanted are never stored. Hutchinson's estimator is unbiased with one line of algebra — and its variance depends on which random vector is used, by a factor that is a property of the matrix, and on a diagonal matrix one choice is exact from the first probe and the other is not.

Eigenvalues, singular values, rank

Keeping the vectors, and losing the bound

Thick restarting keeps the Ritz vectors instead of filtering the starting vector — the same eigenvalues for a third of the products with A. Its residual bound reaches 9.4·10⁻⁴¹ while the residual it bounds sits at 5.7·10⁻⁵, and the eigenvalues are correct to 4.3·10⁻¹⁴ the whole time, so nothing reports it.

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.

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.

Regularisation, and the answer that is chosen

The basis decides what a filter is

The vocabulary of regularisation is spectral — a method keeps a component or discards it, and the weights are a function of the singular value. Row-normalising a symmetric blur so that it preserves a constant makes it 8.6% asymmetric, and that is enough to move GMRES's weights from 7·10⁻¹⁴ off a function of σ to 4.4·10⁻².

Sparsity, and what elimination costs

The factor is not sparse

A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone.

Methods that were designed apart

The step that stops mattering

Regularise the problem the iteration has built rather than the problem it was given, and the error curve stops turning. The unregularised run ends 1,127 times above its own best; the same run with a penalty inside it ends 1.000000000003 times above.

Eigenvalues, singular values, rank

The vector was what was wanted

Nobody who computes a matrix exponential wants the matrix. They want e^{At}b — one vector, the state of a system at a later time. Twenty matrix–vector products get it to sixteen digits on a hundred-by-hundred problem, without ever forming a hundred-by-hundred exponential, and the exponential that does get computed is twenty by twenty.

The whole library · All essays · What must fail