Numerical Linear Algebra
where the algebra stops being the arithmetic
When the index is a tuple

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.
29 min read 25 figures Where the algebra stops being the arithmeticExact ground truthMeasured, not asserted

Worth reading first: Changing the condition number on purpose · An operator with no entries · A block nobody can call sparse.

Three fields of this collection are answers to one sentence: the matrix is too large to store. Sparsity says most of its entries are zero. A matrix-free operator says never form it, and let it be a subroutine. A hierarchical representation says its off-diagonal blocks are numerically low rank, and where they are cut is a decision.

This field is the fourth answer, and it is the only one that changes what an index is.

The 216 eigenvalues of the 3-dimensional Laplacian on 6 points a side, computed against their closed formThe matrix has 216 rows and 46,656 entries, and is a sum of 3 Kronecker products of one 6 × 6 matrix — 108 numbers. Its eigenvalues are every sum of 3 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 7.08·10⁻¹³. The smallest eigenvalue is 0.5942 and the largest 11.41, so the condition number is 19.2 — which is the part the structure does not help with.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
Fig. 1 Every eigenvalue of a matrix with n³ rows, against a closed form that was written down before anything ran. The marks are what a decomposition of the assembled matrix returns; the line is a sum of sines.

The size nobody quotes

A discretisation on a grid of n points along each of d axes has n^d unknowns. That number is quoted constantly and it is not the difficult one. The difficult one is the size of the matrix, which is the square of it: n^{2d}.

At n = 100 and d = 2 that is 10⁸ entries, which is 800 megabytes and merely annoying. At d = 3 it is 10¹², which is eight terabytes. At d = 5 it is 10²⁰, which is not a large number of entries — it is a number of entries with no physical meaning, exceeding the count of atoms in a gram of anything by enough that no improvement in storage will ever matter.

The sparsity field’s answer applies and is the right one for a finite-difference stencil: the matrix has at most 2d + 1 non-zeros a row, so it is 5n² numbers in two dimensions rather than n⁴. That is the answer this collection has used since the two-dimensional model problem arrived, and it is not what this field is about.

What this field is about is a stronger statement, available for the same matrices and for many that are not sparse at all:

The matrix is a sum of d Kronecker products of n × n matrices.

That is d·n² numbers. Nothing has been approximated, no tolerance has been chosen, and no property of the entries has been exploited. The matrix was never anything else. Writing it out as n^{2d} numbers was the mistake, and it is a mistake that is made because the notation for the alternative is not taught alongside the notation for the matrix.

What the product is

For A of size m × n and B of size p × q, the Kronecker product A ⊗ B is the mp × nq matrix built by replacing every entry of A with that entry times the whole of B. It has one row for each pair (i, p) and one column for each pair (j, q), which is the sentence this field turns on: its index is a pair.

The two-dimensional Laplacian on an interior grid is the standard example. Number the unknowns so that u_{i,j} is at position i·n + j, let T = tridiag(−1, 2, −1) be the one-dimensional operator, and the five-point stencil is exactly

A = T ⊗ I + I ⊗ T

The first term differences along the first index and leaves the second alone; the second does the reverse. That is what a Kronecker sum is, and the whole of the two-dimensional operator is in it.

Multiplications in one product with the 2-dimensional Laplacian, assembled against reshapedBoth curves are exact counts rather than estimates. The assembled matrix has n^2 rows, so multiplying by it is n^4 multiplications; applying the 2 one-dimensional factors along their own indices is 2·n^3. The fitted slopes are 4.00 and 3.00 against 4 and 3 exactly. At n = 256 that is 4.29·10⁹ against 3.36·10⁷, a factor of 128 — and the ratio is n^1, so it grows with every size rather than settling.10¹10²10²10⁴10⁶10⁸10¹⁰n, points along one axismultiplicationsassembled: n^4reshaped: 2n^3two exponentsfitted dense slope42d4fitted factored slope3d + 13ratio at n = 256128nothing is approximatedthe matrix was never anything else
Fig. 2 The two-dimensional case, counted. The assembled matrix is n⁴ multiplications per product and the factored one is 2n³, and the ratio is n — so it grows with every refinement of the grid.

For d axes there are d terms and the same construction, with the identity in every position but one. The storage is d·n² whatever d is, which is the property that makes the format worth having: the matrix’s size is exponential in d and its description is linear in it.

The spectrum, which is written down

The one-dimensional Laplacian has its whole spectrum in closed form. Its eigenvalues are

λ_k = 4 sin²(kπ / 2(n + 1)), k = 1 … n

with eigenvectors the discrete sines, and this collection has used them since the iterative field opened, to give the iterative field a rate that is known rather than measured.

The Kronecker sum inherits it exactly. If A has eigenvalues α and B has eigenvalues β, then A ⊗ I + I ⊗ B has eigenvalues α_i + β_j — every sum, one for each pair of indices — and its eigenvectors are the Kronecker products of theirs. So the spectrum of the d-dimensional model problem is every sum of d numbers drawn from that list of sines, and it is known before anything runs.

The hero figure is that statement, checked. The matrix is assembled at a size where a Jacobi decomposition of it can be afforded, its eigenvalues are computed and sorted, and they are compared against the closed form entry by entry. The largest disagreement anywhere is at the rounding level.

The 196 eigenvalues of the 2-dimensional Laplacian on 14 points a side, computed against their closed formThe matrix has 196 rows and 38,416 entries, and is a sum of 2 Kronecker products of one 14 × 14 matrix — 392 numbers. Its eigenvalues are every sum of 2 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 4.01·10⁻¹³. The smallest eigenvalue is 0.08741 and the largest 7.913, so the condition number is 90.52 — which is the part the structure does not help with.033669913216502468eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix196numbers that describe it392λ smallest0.087λ largest7.9worst |computed − exact|4·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 3 The two-dimensional case, where 196 eigenvalues of a 196-row matrix are every sum of two numbers from a list of fourteen.

That is exact-ground-truth in a place it is rarely available. A collection that measures error has to have something to measure against, and its usual recourse is a better approximation — a solve in higher precision, a decomposition of a smaller problem — which is a comparison against a computation rather than against an answer. Here there is an answer, for an operator nobody can store, at every d.

The product, which never needs the matrix

The practical half is a product. Given x with n^d entries, the multiplication (A_1 ⊗ A_2 ⊗ … ⊗ A_d)x is computed by reading x as a d-way array and applying each factor along its own index.

In two dimensions that is one line. Reshape x into an n × n matrix X and

(A ⊗ B) x = vec(A X Bᵀ)

with vec the same row-major reading that turned the grid into a vector in the first place. Two dense matrix multiplications of n × n matrices — 2n³ multiplications — against the assembled route’s n⁴.

The count is the argument, and it is returned by the code beside the answer rather than described. At n = 10 and d = 3 the assembled route is 10⁶ multiplications and the reshaped one is 3·10⁴, and the two answers agree to 4·10⁻¹⁶ relative. At n = 256 and d = 3 the assembled route is 2.8·10¹⁴ and the reshaped one is 1.3·10¹⁰.

Multiplications in one product with the 4-dimensional Laplacian, assembled against reshapedBoth curves are exact counts rather than estimates. The assembled matrix has n^4 rows, so multiplying by it is n^8 multiplications; applying the 4 one-dimensional factors along their own indices is 4·n^5. The fitted slopes are 8.00 and 5.00 against 8 and 5 exactly. At n = 256 that is 1.84·10¹⁹ against 4.4·10¹², a factor of 4.19·10⁶ — and the ratio is n^3, so it grows with every size rather than settling.10¹10²10²10⁶10¹⁰10¹⁴10¹⁸n, points along one axismultiplicationsassembled: n^8reshaped: 4n^5two exponentsfitted dense slope82d8fitted factored slope5d + 15ratio at n = 2564.2·10⁶nothing is approximatedthe matrix was never anything else
Fig. 4 Four indices, where the two exponents are 8 and 5. The gap between the lines is n³ and the picture is the same shape at every d, with the two slopes always d − 1 apart.

The general form of the saving is worth stating because it is the one place in this collection where a difficulty and its remedy scale together. The assembled product costs n^{2d} and the reshaped one costs d·n^{d+1}, so the ratio is n^{d−1}/d. Every index added multiplies the saving by another factor of n. Everywhere else on this site a remedy is a constant factor and the difficulty is not.

What it does not buy

A format that made a difficulty vanish would be worth suspecting, and this one does not.

The condition number is untouched. κ(A ⊗ B) = κ(A)·κ(B) exactly, which is easy to see from the singular values — those of a Kronecker product are all the products of the two sets — and is checked here on a pair of small matrices rather than asserted. Two well-conditioned factors give a product whose condition number is their product, and the d-fold version is a d-th power.

The 81 eigenvalues of the 4-dimensional Laplacian on 3 points a side, computed against their closed formThe matrix has 81 rows and 6,561 entries, and is a sum of 4 Kronecker products of one 3 × 3 matrix — 36 numbers. Its eigenvalues are every sum of 4 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 3.5·10⁻¹³. The smallest eigenvalue is 2.343 and the largest 13.66, so the condition number is 5.828 — which is the part the structure does not help with.01428425670036912eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix81numbers that describe it36λ smallest2.3λ largest14worst |computed − exact|3.5·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 5 Four indices, on a grid of three points a side, which is where the assembled matrix stops being affordable. The condition number of the operator is on the badge and it is not a small number.

For the model problem the arithmetic is worse than that. The one-dimensional Laplacian’s condition number is about (2(n+1)/π)², so it grows like n². The Kronecker sum’s is the ratio of the largest sum to the smallest, and both scale with d — the largest eigenvalue is d times the largest one-dimensional eigenvalue and the smallest is d times the smallest — so the sum has the same condition number as one of its terms, independent of d.

That is a genuinely favourable accident and it is an accident of the sum. A Kronecker product of d copies of T has a condition number that is the d-th power of T’s, and T’s is already 4,135 at n = 100, so the cube of it is 7·10¹⁰. The distinction between the two constructions is not decoration.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 6 What a condition number of that size means for an answer, from the essay that introduced it: the amplification from a backward error to a forward one, and nothing about the storage format changes it.

Where it sits beside the other three

The four answers to the matrix is too large are not alternatives, and the difference between them is worth naming precisely because they are usually all available at once.

Sparsity is a property of the entries. An entry is zero or it is not; the problem decides, and a code either stores the non-zeros or it does not. Nothing is chosen.

A matrix-free operator is a property of the code. The entries may be anything; what is available is the action, and this collection has an essay about what survives when a matrix is only a subroutine. Again nothing is chosen — the problem hands over a subroutine or it does not.

A hierarchical representation is the first one whose cost is chosen. The off-diagonal blocks are numerically low rank, the place their singular values are cut is a number somebody types, and the storage is a function of that number rather than of the matrix.

A Kronecker representation is none of these. It is not an approximation at any accuracy, so there is nothing to choose; it is not a property of the entries, since a Kronecker product of two dense matrices is dense; and it is not a property of the code, since the entries are perfectly available. It is a property of how the problem was written down — of the fact that the domain is a product of intervals and the operator acts along one axis at a time.

CG iterations on the model problem, with the two preconditioners a matrix-free code can and cannot haveAt 256 unknowns, conjugate gradients takes 55 iterations unpreconditioned. Diagonal preconditioning — n probes, or the analytic diagonal, and therefore available without any entries — takes 55, which is the same number, because this operator's diagonal is constant and scaling by it is scaling by a scalar. An incomplete Cholesky takes 21, a factor of 2.62, and it is defined by the sparsity pattern: there is no way to ask a subroutine for it. That factor is the price of an operator with no entries.36861361862360102030405060unknownsCG iterations to 10⁻¹⁰no preconditionerdiagonal (available)incomplete Choleskyat 256 unknownsunpreconditioned55diagonal55incomplete Cholesky21what the pattern buys2.6a Krylov method needs only productsand a preconditioner needs the entries
Fig. 7 The second answer, from the field that measures it: what an operator with no entries costs and what it gives up, which is a different trade from this one entirely.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 8 And the third, from the hierarchy field: storage per unknown that grows like a logarithm because the accuracy was chosen. Nothing on this page has an accuracy to choose.

The failure of the format is a failure of the geometry

The condition under which the representation exists is that the operator separates. A finite-difference Laplacian on a rectangle separates. The same operator on an L-shaped domain does not, because the index set is not a product; a variable coefficient a(x, y) that is not itself a product does not; and an unstructured mesh does not, because there are no axes to speak of.

That is a much sharper restriction than sparsity’s, which is why this field does not replace the others. What it replaces is the class of problems where it applies — tensor-product discretisations, which is most of what a finite-difference or spectral code produces, and separable kernels.

The three unfoldings of a 10 × 11 × 12 wave tensor, and the singular values of eachA tensor has one matrix per index — put that index on the rows and every other index down the columns — and each of those matrices has an ordinary rank. Here they are 2, 2, 2 at a relative tolerance of 10⁻⁸, from a tensor of 1320 entries whose modes are of different lengths. Nothing requires the three numbers to agree, and nothing requires any of them to be the tensor's own rank: they are three different matrices built from one array. The leading singular values are 18.7, 19.3, 18.7, each normalised to its own mode below.02468101210⁻¹⁶10⁻¹²10⁻⁸10⁻⁴1index of the singular valueσ ⁄ σ₁the tolerance the ranks are read atmode 1 · rank 2mode 2 · rank 2mode 3 · rank 2wave: three matrices, one arrayentries1320mode-1 rank2mode-2 rank2mode-3 rank2‖T‖26three ranksand none of them is the tensor's
Fig. 9 A tensor that separates, seen through its three unfoldings: every one of them is rank two, which is what a sum of two products of one-variable functions looks like from any index.
The three unfoldings of a 10 × 11 × 12 noise tensor, and the singular values of eachA tensor has one matrix per index — put that index on the rows and every other index down the columns — and each of those matrices has an ordinary rank. Here they are 10, 11, 12 at a relative tolerance of 10⁻⁸, from a tensor of 1320 entries whose modes are of different lengths. Nothing requires the three numbers to agree, and nothing requires any of them to be the tensor's own rank: they are three different matrices built from one array. The leading singular values are 14.2, 13.6, 13.4, each normalised to its own mode below.02468101210⁻¹⁶10⁻¹²10⁻⁸10⁻⁴1index of the singular valueσ ⁄ σ₁the tolerance the ranks are read atmode 1 · rank 10mode 2 · rank 11mode 3 · rank 12noise: three matrices, one arrayentries1320mode-1 rank10mode-2 rank11mode-3 rank12‖T‖36three ranksand none of them is the tensor's
Fig. 10 And one that does not. The same three matrices, all of full rank, on an array with no structure at all — the counterweight every claim in this field needs.

The refusal

The claim under test is the easiest mistake this field offers, and it is easy because both statements are true of something.

The eigenvalues of a Kronecker product are the products of the eigenvalues. The eigenvalues of a Kronecker sum are the sums. Confusing the two gives a spectrum with the right number of entries, the right sign, and a plausible spread, and it is wrong at every entry.

So the assertion is fed exactly that: the whole spectrum of T ⊗ I + I ⊗ T is computed from the assembled matrix, the set of pairwise products of T’s spectrum is formed, the two sorted lists are compared, and the claim that they agree is required to fail. It does, at the first entry — the smallest sum is 2λ₁ and the smallest product is λ₁², and for the five-point Laplacian at n = 5 those are 0.536 and 0.0719.

The second refusal in the same file is narrower and is about the code rather than the algebra. A mode-k unfolding is a specific matrix, and a routine that folded a mode-1 unfolding back along mode 2 would return an array of exactly the right shape whenever two of the dimensions happen to agree. That is fed a 3 × 3 × 2 array and required to fail.

What the rest of the field is

Everything above is about a matrix. The word tensor has not appeared, because none of it needs one: a Kronecker sum is an ordinary matrix with an index that happens to be a tuple, and its arithmetic is ordinary matrix arithmetic in a different order.

The rest of this field takes the same step one level further, to arrays with three or more indices treated as objects in their own right, and the news there is worse. Every theorem this collection has relied on about rank stops being true: the best low-rank approximation need not exist, the rank depends on the field the entries are read over, and there is no decomposition that is orthogonal and diagonal at once. What survives is not the definition of a decomposition but the algorithm — take the singular value decompositions of the reshapes — and existence, computability and a quasi-optimality factor all come back with it.

A rank-two sequence approaching a rank-three tensor: the residual falls like 1/n and the terms grow like nA_n = n(e₁ + e₂/n)⊗³ − n·e₁⊗³ has rank two for every n and converges to a tensor of rank three. The falling curve is ‖A_n − A‖, which is √(3/n² + 1/n⁴) exactly and reaches 0.00677 at n = 256; the rising one is the norm of the larger of its two rank-one terms, 256. Their product runs 2.519, 1.917, 1.777 … 1.73209, descending onto √3 = 1.73205. So getting one digit closer costs a factor of ten in the size of the pieces, for ever, and the infimum of the distance is zero while no rank-two tensor attains it.110¹10²10⁻²10⁻¹110¹10²n‖A_n − A‖ and the largest term's norm‖A_n − A‖the larger of its two termsan infimum that is not attainedn256‖A_n − A‖0.0068largest term256their product1.7√31.7the distance goes to zeroand nothing reaches it
Fig. 11 Where that goes, from two essays on: a sequence of rank-two tensors converging to a tensor of rank three, with nothing at the limit for it to converge to.
The truncation error of a smooth tensor against the rank kept, between the two bounds the theorem givesThe middle curve is the measured error of the projection; the upper dashed one is √(Σ_k tail_k²), which the theorem says it cannot exceed, and the lower one is max_k tail_k, which the best possible error cannot fall below. They are a factor of √3 apart. The measurement is that the projection sits on the upper one, and not between them: the ratio of error to bound runs 0.688, 0.717, 0.879, 0.933 … 0.999982, so by rank 10 the bound is attained to five decimals and the ratio to the lower bound is 1.7320 against √3 = 1.7321. That reads as a bad result and is not one — what it says is that the lower bound is weak, which only a second measurement can establish.024681010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnssmooth: pinned to the upper boundrank 10 error1.1·10⁻¹¹its upper bound1.1·10⁻¹¹the lower bound6.3·10⁻¹²error ⁄ bound1error ⁄ lower1.7inside the boundand sitting on it
Fig. 12 And how it is repaired: the projection whose error sits between two computable bounds a factor of √d apart, and on the upper one.

At other settings

The 216 eigenvalues of the 1-dimensional Laplacian on 216 points a side, computed against their closed formThe matrix has 216 rows and 46,656 entries, and is a sum of 1 Kronecker products of one 216 × 216 matrix — 46656 numbers. Its eigenvalues are every sum of 1 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 2.62·10⁻¹³. The smallest eigenvalue is 0.0002096 and the largest 4, so the condition number is 19080 — which is the part the structure does not help with.0367210814418021601234eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix216numbers that describe it4.7·10⁴λ smallest2.1·10⁻⁴λ largest4worst |computed − exact|2.6·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 13 One index, where the Kronecker sum is the one-dimensional operator itself and the closed form is the sine formula unmodified — the base case the other three stops are built from.
Multiplications in one product with the 3-dimensional Laplacian, assembled against reshapedBoth curves are exact counts rather than estimates. The assembled matrix has n^3 rows, so multiplying by it is n^6 multiplications; applying the 3 one-dimensional factors along their own indices is 3·n^4. The fitted slopes are 6.00 and 4.00 against 6 and 4 exactly. At n = 256 that is 2.81·10¹⁴ against 1.29·10¹⁰, a factor of 2.18·10⁴ — and the ratio is n^2, so it grows with every size rather than settling.10¹10²10²10⁵10⁸10¹¹10¹⁴n, points along one axismultiplicationsassembled: n^6reshaped: 3n^4two exponentsfitted dense slope62d6fitted factored slope4d + 14ratio at n = 2562.2·10⁴nothing is approximatedthe matrix was never anything else
Fig. 14 Three indices: exponents of 6 and 4, and a ratio at n = 256 of 4.3·10⁴.
Multiplications in one product with the 6-dimensional Laplacian, assembled against reshapedBoth curves are exact counts rather than estimates. The assembled matrix has n^6 rows, so multiplying by it is n^12 multiplications; applying the 6 one-dimensional factors along their own indices is 6·n^7. The fitted slopes are 12.00 and 7.00 against 12 and 7 exactly. At n = 256 that is 7.92·10²⁸ against 4.32·10¹⁷, a factor of 1.83·10¹¹ — and the ratio is n^5, so it grows with every size rather than settling.10¹10²10²10⁷10¹²10¹⁷10²²10²⁷n, points along one axismultiplicationsassembled: n^12reshaped: 6n^7two exponentsfitted dense slope122d12fitted factored slope7d + 17ratio at n = 2561.8·10¹¹nothing is approximatedthe matrix was never anything else
Fig. 15 Six, where the assembled count reaches 10²⁸ at the right-hand end of the axis and is drawn only so that the other line has something to be compared against.
The three unfoldings of a 10 × 11 × 12 smooth tensor, and the singular values of eachA tensor has one matrix per index — put that index on the rows and every other index down the columns — and each of those matrices has an ordinary rank. Here they are 8, 8, 8 at a relative tolerance of 10⁻⁸, from a tensor of 1320 entries whose modes are of different lengths. Nothing requires the three numbers to agree, and nothing requires any of them to be the tensor's own rank: they are three different matrices built from one array. The leading singular values are 2.37, 2.37, 2.37, each normalised to its own mode below.02468101210⁻¹⁶10⁻¹²10⁻⁸10⁻⁴1index of the singular valueσ ⁄ σ₁the tolerance the ranks are read atmode 1 · rank 8mode 2 · rank 8mode 3 · rank 8smooth: three matrices, one arrayentries1320mode-1 rank8mode-2 rank8mode-3 rank8‖T‖2.4three ranksand none of them is the tensor's
Fig. 16 The three unfoldings of a smooth array of different mode lengths, where the ranks are close and are not equal — three matrices built from one array, with no reason to agree.
The three unfoldings of a 10 × 11 × 12 hilbert tensor, and the singular values of eachA tensor has one matrix per index — put that index on the rows and every other index down the columns — and each of those matrices has an ordinary rank. Here they are 8, 8, 8 at a relative tolerance of 10⁻⁸, from a tensor of 1320 entries whose modes are of different lengths. Nothing requires the three numbers to agree, and nothing requires any of them to be the tensor's own rank: they are three different matrices built from one array. The leading singular values are 3.2, 3.2, 3.2, each normalised to its own mode below.02468101210⁻¹⁶10⁻¹²10⁻⁸10⁻⁴1index of the singular valueσ ⁄ σ₁the tolerance the ranks are read atmode 1 · rank 8mode 2 · rank 8mode 3 · rank 8hilbert: three matrices, one arrayentries1320mode-1 rank8mode-2 rank8mode-3 rank8‖T‖3.3three ranksand none of them is the tensor's
Fig. 17 And the d-way version of the matrix this site has used for exact ground truth from the beginning.
Numbers stored against the number of indices, at n = 20 and rank 4: the tensor, its core, and a trainThe tensor is n^d, which at d = 20 is 1.05·10²⁶. A Tucker representation of it is r^d + d·n·r — the core is still exponential in d, so at rank 4 it is 1.1·10¹², smaller than the tensor by 9.54·10¹³ and still unstorable. A train is (d − 2)nr² + 2nr, which is 5,920 — linear in d. The three lines are the field's whole argument: fixing the definition of the decomposition does not fix the size of what it returns, and the second fix is the same projection cut in a different place.15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core1.1·10¹²train5920core ⁄ train1.9·10⁸entries ⁄ core9.5·10¹³the definition is repairedthe size is not
Fig. 18 What the field’s other formats cost in the same units, for reading against the d·n² above.
Entries against numbers stored, for sin of a sum on 6 points a side, as indices are addedThe upper line is the tensor: 6^d entries, which is a straight line on a logarithmic axis and reaches 46,656 at d = 6. The lower one is the train, which for this family is 4n(d − 1) exactly — 24, 48, 72, 96, 120 — a straight line on a *linear* axis and therefore a logarithm on this one. Its fitted slope against d is 24.0, which is 4n. The two are the same object to within 1.39·10⁻¹⁵, so nothing has been given up: the ratio at d = 6 is 389, and it grows by a factor of n with every index added.123456710¹10²10³10⁴10⁵number of indicesnumbersentries: 6^dstored: 4n(d − 1)exponential against linearentries at d = 64.7·10⁴numbers stored120ratio389slope against d24‖T − T_tt‖ ⁄ ‖T‖1.4·10⁻¹⁵one line is n^dthe other is a constant per index
Fig. 19 And the format that is linear in the number of indices, which is where this field ends up.
The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 20 The first answer, from the sparsity field: the same operator stored as its non-zeros, and what an elimination does to that pattern — which is the comparison this whole page is against.
What one weighted Jacobi sweep does to each mode, at ω = 0.67A curve of the damping factor against mode frequency, falling from one at the left through zero and towards minus one at the right. The right-hand half of the plot is shaded, and two horizontal lines mark the worst damping inside it.00.250.50.751-1-0.500.51mode frequency θ / πdamping factor per sweeppredicted±0.334measureda coarse grid seestwo routes to one factorsmoothing factor, scanned0.33smoothing factor, closed form0.33worst mode disagreement4.4·10⁻¹⁶63 interior points, one sweepthe left-hand end is what the coarse grid is for
Fig. 21 The eigenvectors the closed form comes with, from the multigrid field: discrete sines, whose Kronecker products are the eigenvectors of everything above, and what a smoother does to each of them.
A 4-level hierarchy on 15 points, and the coarse operator it impliesRows of dots, each row half the length of the one above it, with lines joining every coarse point to the three fine points it interpolates to.15731pointsfinest gridone unknown — the recursion bottoms out in a divisionthe coarse operator, two ways‖RA_hP − A_2h‖/‖A_2h‖10⁻¹⁸unknowns / finest grid1.7cycle cost, in fine sweeps12each coarse point reaches three fine ones½, 1, ½ — and the restriction is its transpose
Fig. 22 The hierarchy the same operator supports, for scale: four grids of a problem whose matrix on the finest of them is the object this page has been describing.
Jacobi, Gauss–Seidel and SOR at ω = 1.777A semi-logarithmic plot of relative residual against iteration for three stationary methods, with dashed reference curves showing the rate each is predicted to contract at.015030045060075090010⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖JacobiGauss–SeidelSOR ω=1.78closed form vs measuredρ Jacobi, exact0.99ρ measured0.99ρ Gauss–Seidel, exact0.981D Laplacian, n = 24ω optimal = 1.777
Fig. 23 And the rates the closed-form spectrum gives in advance, which is what having an exact spectrum is worth to the field that uses it most.
The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix17 entriestip eliminated first45 entriestip eliminated last17 entries‖A − LLᵀ‖/‖A‖, tip first1.1·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last1.1·10⁻¹⁶dense factor is n(n+1)/2 = 45 · sparse factor is 2n − 1 = 17one row swapped to the endnothing numerical chose between them
Fig. 24 The sparsity field’s other picture of an index order deciding a cost, for reading against the reshape above: the same matrix numbered two ways.
Convergence factor against grid size, with 2 pre-smoothing sweepsThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.003Jacobi at n = 1270.99work exponent, fitted0.079the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 25 And what the field that uses this operator most does with it, for scale: a rate that does not depend on the size, measured on the matrix this page has been describing.

Where this goes next

10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice 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. the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding 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. 051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns A block nobody can call sparse A 96 × 96 block of a kernel matrix has ninety-six nonzero singular values and five that matter. It has no zero entries, it is not described by fewer numbers than it contains, and neither of the two ways this collection already knows to make a large matrix affordable applies to it.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Condition numberCurse of dimensionalityDiscrete laplacianExact ground truthKronecker productKronecker sumMatrix-freeModel problemSeparabilityUnfolding