Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

An iterate that must be made smaller

Applying a Kronecker-sum operator to a low-rank iterate multiplies its ranks by d and adding two of them adds their ranks, so a solver in a compressed format cannot keep what it produces. Every step is followed by a truncation — and whether that truncation is a floor on the residual depends on the right-hand side rather than on the truncation.
33 min read 27 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: The rate the condition number predicts · The format that does not notice the dimension · The rounding that was not the problem.

Every solver in this field produces a sequence of vectors. When the vector is stored densely that is a detail; when it is stored in a compressed format it is the whole problem, because the format is not closed under the operations the method performs.

Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 3The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 3.85·10⁻⁴. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 3. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 3solid: the budgeted residual · light: no budgetbudget 3rank asked for5rank kept3residual, budgeted3.8·10⁻⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored150the step asks for moreat every step
Fig. 1 A conjugate gradient run whose every iterate is cut back to a rank budget. The two step curves near the top are what the method asks for and what it is allowed to keep.

What the arithmetic does to the ranks

Two operations, and both of them grow the representation.

Applying the operator. The d-dimensional model problem is a Kronecker sum, so applying it is d mode products summed. Each mode product leaves the ranks alone and the sum of d of them multiplies them by at most d. So Ax has ranks up to d times x’s.

Adding two iterates. A conjugate gradient step is x ← x + αp and p ← r + βp, and adding two representations concatenates their factors: the ranks add.

Put the two together and one step takes ranks r to something like (d + 1)r. After five steps that is 1,000r on the three-dimensional problem, and the representation is larger than the tensor it is representing.

So a solver in this format has no choice: every step is followed by a truncation back to something affordable, and the whole essay is about what that costs.

Train ranks at each of the 4 cuts of a 5-index tensor on 6 points a side, four familiesCut the index list after position k, put the first k indices on the rows and the rest on the columns, and take the rank of the matrix that results. There are 4 such cuts and each rank is an ordinary matrix rank. sin(x₁ + … + x_d) has rank exactly two at every one of them, for every d, because the addition formula separates it into two terms at every cut — a rank written down rather than measured, and the computed values are 2, 2, 2, 2. The reciprocal family climbs to 9, the product family is one everywhere, and independent normal entries reach 36, which is the largest rank the cut allows. Storage runs sinsum 96, reciprocal 1,206, product 30, noise 10,440 against 7,776 entries.012345061218243036cut after index krank of the reshapesinsum · 2 2 2 2reciprocal · 6 9 9 6product · 1 1 1 1noise · 6 36 36 64 cuts, 4 rankssinsum stored96reciprocal stored1206product stored30noise stored10⁴entries7776one rank per cutand one of them is a theorem
Fig. 2 The quantity being budgeted, from the field that introduces it: one rank per place the index list can be cut.

The measurement the format cannot make

The interesting number is what the step wanted, and a code in this format has thrown it away by the time anyone could ask.

So the arithmetic here is dense, deliberately. The iterate is carried as a full array, its true cut ranks are read at every step, and the truncation is applied by a round trip through the compressed format. That is not how a real solver works and it is the only way to see the quantity the truncation exists to control.

On the three-dimensional problem at ten points a side with a constant right-hand side, the un-truncated step asks for five at every step from the third onwards. A budget of three gives three. The truncation is therefore not an occasional tidy-up; it is happening at every step of every run, and the distance between the two residual curves on the hero figure is what it costs.

Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 1The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 0.119. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 1. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 1solid: the budgeted residual · light: no budgetbudget 1rank asked for5rank kept1residual, budgeted0.12residual, unbudgeted1.6·10⁻¹⁴numbers stored30the step asks for moreat every step
Fig. 3 The tightest budget drawn, where the iterate is held at rank one and the residual stalls two decades above where the method would otherwise reach.
Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 5The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 1.61·10⁻¹⁴. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 5. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 5solid: the budgeted residual · light: no budgetbudget 5rank asked for5rank kept5residual, budgeted1.6·10⁻¹⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored350the step asks for moreat every step
Fig. 4 And the budget at which the two residual curves become indistinguishable, which is not a property of the method or the tolerance.

What the budget is worth

The ladder is the measurement, and it has to be run twice because the answer depends on something other than the truncation.

With a constant right-hand side, the floors are

0.118, 0.0107, 3.6·10⁻⁴, 5.4·10⁻⁶, 1.1·10⁻¹³

at budgets of one, two, three, four and six. About two decades a column, and then the floor disappears — because the solution of this problem is a train of rank five, and once the budget reaches five the truncation has stopped removing anything.

With independent normal entries in the right-hand side, the same ladder on the same operator gives

0.969, 0.925, 0.874, 0.801, 0.656, 0.470

at budgets of one to eight out of a possible ten. Eight columns of a ten-column answer buy a factor of two.

So a floor is not a property of the truncation. It is the distance from the answer to the set the truncation projects onto, and the same budget is free on one problem and useless on the next.

The residual a rank budget leaves, against the budget, on two right-hand sides of the same operatorBoth runs use the same method on the same matrix for 70 steps and differ only in b. With a constant right-hand side the floor falls 0.118, 0.0107, 3.61·10⁻⁴, 5.4·10⁻⁶, 1.07·10⁻¹³, 1.07·10⁻¹³, 1.07·10⁻¹³ — about two decades per column — and disappears at rank 5, because the solution *is* a train of that rank and the truncation has stopped removing anything. With independent normal entries the same ladder gives 0.969, 0.925, 0.874, 0.801, 0.735, 0.656, 0.47: 8 columns of a 10-column answer buy a factor of 2.06. So a floor is not a property of the truncation. It is the distance from the answer to the set the truncation projects onto, and the same budget is free on one problem and useless on the next.012345678910⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹rank budgetresidual the run stalls atb with no structure at allb constant: the answer is a traintwo ladders, one truncationstructured, rank 10.12structured, rank 45.4·10⁻⁶structured, rank 81.1·10⁻¹³unstructured, rank 10.97unstructured, rank 80.47the floor is not the truncation'sit is the answer's
Fig. 5 The two ladders on one pair of axes. Same operator, same method, same budgets — and the only difference is what is being solved for.

The step, in detail

It is worth writing out what one step actually does to the representation, because the growth is not a vague statement about complexity.

Start with x of rank r. Then:

  • Ap has rank at most d·r_p, because the operator is a sum of d terms each of which preserves ranks.
  • αp has rank r_p, since scaling does not change a factorisation.
  • x + αp has rank at most r + r_p.
  • b − Ax has rank at most r_b + d·r.
  • r + βp has rank at most (r_b + d·r) + r_p.

Every one of those is an upper bound and every one of them is attained generically, which the measured before column on the hero figure confirms: the un-truncated cut ranks climb to the ceiling the reshape itself imposes within three steps.

The truncation applied afterwards is the format’s own, so it costs d − 1 decompositions of reshapes of the iterate. On a problem where the iterate is genuinely low rank those decompositions are of small matrices and the cost is negligible; on a problem where it is not, the truncation is the dominant expense and it is buying nothing, which is the second ladder above.

The singular values of a sum of two rank-4 blocks, and the 4 a truncation has to discardA rank-4 block times a vector is a rank-4 block times a vector. A rank-4 block times a rank-4 block is a rank-4 block. A rank-4 block PLUS a rank-4 block is a rank-8 block, exactly, and the 8 bars here are why: the sum of two 4-dimensional spaces is generally 8-dimensional, and none of the 8 singular values is small. Truncating back to 4 costs 61.9 per cent of the block. Below the 8th the values are the unit roundoff, which is the check that the doubling is exact rather than approximate. Every product, every factorisation and every Schur complement inside this format is a chain of these, and there is nothing else to do: without the truncation the ranks double at every level and the format is dense by the bottom.σ ⁄ σ₁ of the sum, 64 × 64σ1, kept1σ2, kept0.922σ3, kept0.812σ4, kept0.785σ5, discarded0.778σ6, discarded0.758σ7, discarded0.658σ8, discarded0.571σ92.63·10⁻¹⁶σ102.3·10⁻¹⁶the operation that is not closedrank of each term4rank of the sum8truncated back to4cost of the truncation0.62the best there is0.62two planesmake a four-space
Fig. 6 The hierarchy field’s version of the same arithmetic, from the essay that measures it: a rank that doubles under addition, exactly.

Why that is not what the neighbouring fields found

This collection has two other essays about a truncation inside an iteration, and both of them find a floor. Putting the three together is the point of this one.

The hierarchy field performs a Cholesky entirely inside a compressed format, with the number of truncations a count rather than a parameter — ninety-eight of them at the deepest level — and finds the residual below the representation’s own error at every depth. Its conclusion is that the leaf size can be chosen from the machine, because there is no accuracy term in the trade.

The residual-gap essay finds a conjugate gradient run whose reported residual falls below the unit roundoff while its true residual does not, and traces the gap to the size of the iterates rather than to any single rounding.

Both of those are about a truncation whose error accumulates or fails to. This page is about a truncation whose error may be zero — and the condition for that is not a property of the truncation at all. That is a genuinely different finding, and it is the one a practitioner needs, because it says the question to ask about a rank budget is not how much accuracy does it cost but is the answer in the set.

‖A − LLᵀ‖ ⁄ ‖A‖ for a Cholesky computed in the format, against how many truncations it performedThe recursion is the textbook one and every step stays in the format: the off-diagonal factor is exact, because a triangular solve against one factor of a rank-k block leaves a rank-k block, and the Schur complement is a rank-k addition to a hierarchical matrix, which is the only approximate step there is. So the number of truncations is a count — one per admissible block in the subtree being updated, at every level — and it runs from 0 at a leaf of 128 to 98 at a leaf of 8. The residual goes from 3.856·10⁻¹² to 1.035·10⁻¹¹, and all of that movement is the representation, whose own error is the upper curve and moves by the same factor: the ratio between them is 1.00, 0.93, 0.93, 0.91, 0.86, never above one. A hundred approximate operations contributed nothing measurable to the answer.02040608010010⁻¹²10⁻¹¹10⁻¹⁰truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations3.9·10⁻¹²residual, 98 truncations10⁻¹¹representation, deepest1.2·10⁻¹¹residual ⁄ representation0.86levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 7 The hierarchy field’s version, from the essay that measures it: ninety-eight truncations and a residual that does not notice them.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹² and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 1.6·10¹⁰ apart, with the reported one at 4.01·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01326395265789110410⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best4·10⁻²¹the answer's, at its best6.2·10⁻¹¹unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹¹iterations drawn105the recurrence remembers every roundingand the stopping test is written in it
Fig. 8 And the residual-gap essay’s, where the quantity a method reports and the quantity it has separate.

What the orthogonality costs

The method here is conjugate gradients, and truncating inside it breaks the thing the method is built on.

A short recurrence works because each new direction is automatically orthogonal to every previous one, which is a consequence of exact arithmetic and of the iterates lying in the Krylov subspace. Projecting each iterate onto a low-rank set takes it out of that subspace, so the conjugacy is lost immediately and completely — not gradually, as it is in finite precision.

Nothing here repairs it. The residual is recomputed from the iterate at every step rather than updated by the recurrence, so the number plotted is a true residual of a real vector, and the method’s own bookkeeping is left to be as wrong as the truncation makes it. That is the honest arrangement: a plot of the recurrence’s own residual would show something much better than the truth, which is exactly the defect the residual-gap essay is about.

What a real code does about it is restart, or use a method that does not rely on a short recurrence, or re-orthogonalise against a stored basis — and all three cost storage, which is what the format was for.

The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 9 The property being broken, from the essay that measures it: an orthogonality that is a number rather than an assumption.
Conjugate gradients at κ = 103, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound85 steps40×40, spectrum spread evenly in logbound permits 448
Fig. 10 And the rate the method has when nothing is truncated, which is what the budgeted curves are measured against.

What decides whether a right-hand side has a rank

The two ladders above are extremes and the practical question is which one a given problem resembles.

For the model problem the answer is available in closed form and is worth having. The solution is A⁻¹b, and the inverse of a Kronecker sum is not a Kronecker sum but is within any accuracy of a short sum of them — half a term a decade, which the tensor field measures. So A⁻¹ applied to a rank-one b gives something of low rank, and applied to a full-rank b gives something of full rank.

That is a statement about b and not about A. A right-hand side that is a product of one-variable functions, or a sum of a few such products, gives a low-rank solution; a right-hand side sampled from data gives whatever it gives. And the rank of b is cheap to compute before any solving happens: it is d − 1 decompositions of reshapes.

So the honest workflow is to measure the right-hand side’s ranks first, and to treat a large one as evidence that the format is the wrong one rather than as a reason to raise the budget.

The Kronecker spectrum of the inverse of a Kronecker sum on 8 points a sideA Kronecker product B ⊗ C, read as a four-index array and cut between its two index pairs, is exactly rank one. The inverse of T ⊕ T is rank 8 at the same cut, so the format is not closed under inversion — which is why a solve in it goes through the eigenbasis rather than through an inverse. What the spectrum says is that it is nearly closed: the singular values are 1, 0.19, 0.0262, 0.00214 of the first, and the number of Kronecker terms needed runs 3, 5, 6, 7, 8, 8 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is 0.50 terms a decade — the same shape, and nearly the same number, as the 0.554 columns a decade the hierarchy field measures for a kernel block, arrived at from a different direction entirely.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
Fig. 11 The closed-form half of that, from the tensor field: the inverse of a Kronecker sum, and the number of terms an accuracy buys.
The largest train rank of a 5-index reciprocal tensor against the accuracy asked forEvery point is a decomposition that met its tolerance: the measured errors are 0.0041, 4.2·10⁻⁴, 3.4·10⁻⁵, 3.3·10⁻⁷, 5.7·10⁻⁹, 6.6·10⁻¹², 1.6·10⁻¹³ against tolerances of 10⁻², 10⁻³, 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². The rank runs 3, 4, 5, 7, 8, 10, 11 and the storage 264, 448, 680, 1,160, 1,664, 2,208, 2,504, against 32,768 entries. That is 0.80 of rank per decade of accuracy over the whole range — a constant, with no cliff and no regime where a digit costs more than the last one, which is the same shape the hierarchy field measured for a kernel matrix and is not something either field's geometry promises.-13-11-9-7-5-3-1024681012log₁₀ of the accuracy asked forlargest train rank0.80 of rank per decadea cost that is typed inentries3.3·10⁴stored at 10⁻²264stored at 10⁻¹²2504rank per decade0.8worst error ⁄ tolerance0.57the storage is chosena constant of rank a decade
Fig. 12 And the cost of measuring a right-hand side’s ranks, which is the same sweep at a different tolerance.

Where this leaves the field’s other methods

Three of them transfer and one does not.

Stationary methods transfer cleanly. A Jacobi or Gauss–Seidel sweep on a Kronecker sum is a mode product and a scaling, its rank growth is bounded per step, and there is no recurrence to break. The rates are the ones the model problem’s closed-form spectrum gives, unchanged.

Multigrid transfers, and the coarse grids are cheaper in this format than in any other, because restriction and prolongation are mode products. What does not transfer is the argument that the work per unknown is constant, since the ranks may grow as the grid refines.

Preconditioning transfers and is where the format is strongest: a Kronecker-sum preconditioner is applied by the exact solve the tensor field describes, and applying it does not grow ranks at all.

Krylov methods transfer worst, for the reason two sections above. The method whose whole advantage is a short recurrence is the method a projection breaks.

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. 13 The rates that transfer unchanged, from the field’s own essay: known in advance, on the operator this page solves.
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. 14 And the method whose transfer is partial, whose rate does not notice the size until the ranks do.
Incomplete Cholesky on the 15×15 grid: κ 103 → 9.96A semi-logarithmic plot of relative residual against iteration for conjugate gradients with and without an incomplete Cholesky preconditioner, the preconditioned curve falling faster.071421283542495610⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)103κ(L⁻¹AL⁻ᵀ)10‖A − LLᵀ‖/‖A‖0.0862D Laplacian, n = 225√κ ratio predicts 3.22×
Fig. 15 The one that is strongest here, measured in the units this field prices everything in.

What it would take to do better

Three things a real code does that this measurement does not, and each of them is a different trade.

Truncate less often. The growth is bounded per step, so a code can afford to let the ranks rise for a few steps and truncate every third one. That reduces the number of decompositions and increases the peak storage, and where the optimum sits is a property of the machine rather than of the mathematics — which is the same shape of answer the cost field gives about block sizes.

Truncate adaptively. A fixed budget is the crudest rule. Cutting each iterate at a tolerance relative to the current residual keeps the truncation error below the iteration error at every step, which is the inexact-Newton argument the sequence field makes about inner tolerances: solve to the accuracy the outer loop can use, and no further.

Restart. Since the conjugacy is lost anyway, a method that does not depend on it — a restarted minimal-residual iteration, or plain steepest descent with an optimal step — loses nothing by the truncation and is easier to reason about. The measurement here uses conjugate gradients precisely because it is the case where the damage is clearest.

None of the three changes the finding. All of them are about how efficiently the budget is spent, and the finding is about whether spending it buys anything at all.

A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.507 from a root that is known exactly by construction. The same linear system is solved to relative residuals from 0.3 down to 10⁻¹⁴, at 78 and 1005 conjugate gradient iterations, and the resulting point is 0.2069 and 0.1787 from the root. The curve is flat below about 10⁻³: the linearisation is wrong at second order, so the step cannot land closer than the square of the distance it started at — 0.2572 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻²10⁻¹110¹inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.257distance before the step, 0.507782576211005the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.51its square0.26where η = 10⁻³ lands0.18where η = 10⁻¹⁴ lands0.18iterations for the first257iterations for the second1005the accuracy that is thrown awaymeasured against a root that is known
Fig. 16 The second option’s argument, from the sequence field: an inner accuracy read off the outer residual rather than fixed in advance.
Words moved against the block size, n = 48, M = 100A U-shaped curve on logarithmic axes. Too small a block refactors the panel too often; too large a block does not fit in the fast memory and the tiled update evicts what it is about to read. The minimum is at b = 8, which is √M − 2. The three-line count says √(M/3) = 5, which has the right scaling and the wrong constant.110¹10⁴10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.6·10⁴at b = 13.9·10⁴at b = 243.9·10⁴derived from M with no measurement, and scannedthe two agree
Fig. 17 And the first option’s, from the cost field: an optimum that is a property of the machine.

The counting that makes the format worth it anyway

It is worth ending on the arithmetic that motivates the whole exercise, because two negative measurements in a row misrepresent the situation.

At d = 5 and n = 20 the problem has 3.2 million unknowns, which is a dense vector of 26 megabytes and a matrix that cannot be written down. A train of rank ten stores 5·20·100 = 10,000 numbers — 80 kilobytes — and every operation on it is a handful of small dense multiplications.

The measurement on this page says that the 80 kilobytes buy the answer when the answer is in the format and buy a factor of two when it is not. Both halves are worth having, and the second is worth having because nothing else reports it: a code that ran the second ladder and looked only at its residual would conclude that the problem is hard, when what has happened is that the representation is wrong for it.

The check is cheap and it is available before any solving happens. Decompose the right-hand side, read its cut ranks, and compare them against the budget. A right-hand side that is already at the ceiling has a solution that will be, and no amount of iteration changes that.

Train ranks at each of the 3 cuts of a 4-index tensor on 6 points a side, four familiesCut the index list after position k, put the first k indices on the rows and the rest on the columns, and take the rank of the matrix that results. There are 3 such cuts and each rank is an ordinary matrix rank. sin(x₁ + … + x_d) has rank exactly two at every one of them, for every d, because the addition formula separates it into two terms at every cut — a rank written down rather than measured, and the computed values are 2, 2, 2. The reciprocal family climbs to 9, the product family is one everywhere, and independent normal entries reach 36, which is the largest rank the cut allows. Storage runs sinsum 72, reciprocal 720, product 24, noise 2,664 against 1,296 entries.01234061218243036cut after index krank of the reshapesinsum · 2 2 2reciprocal · 6 9 6product · 1 1 1noise · 6 36 63 cuts, 3 rankssinsum stored72reciprocal stored720product stored24noise stored2664entries1296one rank per cutand one of them is a theorem
Fig. 18 The check, which is d − 1 decompositions of reshapes and costs a fraction of one step of the solver.

The refusal

The claim under test is the one a rank budget invites: that it is a cost rather than an accuracy — spend more memory, get closer.

The assertion that a budget of eight gets the residual below 10⁻⁶ is fed the random right-hand side, on an operator whose answer has ten columns. It stalls at 0.47 and the assertion fails.

That is the right direction for the refusal. Asserting that the budget does floor the residual would pass on the structured right-hand side and on the unstructured one and would say nothing; asserting that it reaches the answer is the flattering reading, and feeding it the case where it is false is what makes the two ladders a finding rather than an illustration.

The file’s other refusals cover the neighbours. One is fed the sin family and required to refuse the claim that a function of a sum is separable into one term. The other is fed a tensor of independent normal entries and required to refuse the claim that the format compresses it — the counterweight this whole area needs.

At other settings

Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 2The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 0.0112. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 2. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 2solid: the budgeted residual · light: no budgetbudget 2rank asked for5rank kept2residual, budgeted0.011residual, unbudgeted1.6·10⁻¹⁴numbers stored80the step asks for moreat every step
Fig. 19 A budget of two, between the two settings the argument is made at.
Conjugate gradients on the three-dimensional model problem with every iterate cut to rank 8The falling curves are relative residuals: the lower one is the same method with no budget, which reaches 1.61·10⁻¹⁴ in 40 steps, and the upper one is the budgeted run, which stops at 1.61·10⁻¹⁴. The two step curves near the top are ranks, on their own scale: the un-truncated step asks for 5 at every step from the third onwards and the budget allows 8. The truncation is therefore not an occasional tidy-up — it is happening at every step, and the distance between the two residual curves is what it costs. The solution of this problem is itself a train of rank five, so a budget of five or more removes nothing and the two curves coincide.081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 8solid: the budgeted residual · light: no budgetbudget 8rank asked for5rank kept8residual, budgeted1.6·10⁻¹⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored350the step asks for moreat every step
Fig. 20 And a budget above the answer’s own rank, where the truncation removes nothing and costs nothing.
Train ranks at each of the 2 cuts of a 3-index tensor on 6 points a side, four familiesCut the index list after position k, put the first k indices on the rows and the rest on the columns, and take the rank of the matrix that results. There are 2 such cuts and each rank is an ordinary matrix rank. sin(x₁ + … + x_d) has rank exactly two at every one of them, for every d, because the addition formula separates it into two terms at every cut — a rank written down rather than measured, and the computed values are 2, 2. The reciprocal family climbs to 6, the product family is one everywhere, and independent normal entries reach 6, which is the largest rank the cut allows. Storage runs sinsum 48, reciprocal 288, product 18, noise 288 against 216 entries.012301234567cut after index krank of the reshapesinsum · 2 2reciprocal · 6 6product · 1 1noise · 6 62 cuts, 2 rankssinsum stored48reciprocal stored288product stored18noise stored288entries216one rank per cutand one of them is a theorem
Fig. 21 The ranks of the arrays involved, at the number of indices this page solves at.
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. 22 And what the format stores, which is the reason to put up with any of this.
Multiplications in a 3-dimensional model solve, through the eigenbasis against a dense factorisationA Kronecker sum's eigenvectors are the Kronecker products of its factors' eigenvectors, so the change of basis that diagonalises an operator with 1,728 rows is 3 changes of basis along 3 indices. The whole solve is a transform, 1,728 divisions and a transform back: 1.24·10⁵ multiplications at n = 12, against a dense factorisation's 3.44·10⁹, a factor of 2.76·10⁴. The fitted exponent is 4.00 against 3d = 9. The only decompositions taken are of the 3 one-dimensional factors, which on the model problem are the same matrix — so there is exactly one, of size 12. Every solve reproduces its right-hand side to 8.36·10⁻¹⁵.10¹10²10⁴10⁶10⁸10¹⁰n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis3 decompositions of an n × nunknowns1728multiplications1.2·10⁵dense factorisation3.4·10⁹fitted exponent4‖Ax − b‖ ⁄ ‖b‖8.4·10⁻¹⁵nothing of size n^dis ever factorised
Fig. 23 The direct alternative on the same operator, from the tensor field: a solve through the eigenbasis, which needs no iteration and no truncation at all.
The singular values of a sum of two rank-4 blocks, and the 4 a truncation has to discardA rank-4 block times a vector is a rank-4 block times a vector. A rank-4 block times a rank-4 block is a rank-4 block. A rank-4 block PLUS a rank-4 block is a rank-8 block, exactly, and the 8 bars here are why: the sum of two 4-dimensional spaces is generally 8-dimensional, and none of the 8 singular values is small. Truncating back to 4 costs 61.9 per cent of the block. Below the 8th the values are the unit roundoff, which is the check that the doubling is exact rather than approximate. Every product, every factorisation and every Schur complement inside this format is a chain of these, and there is nothing else to do: without the truncation the ranks double at every level and the format is dense by the bottom.σ ⁄ σ₁ of the sum, 64 × 64σ1, kept1σ2, kept0.922σ3, kept0.812σ4, kept0.785σ5, discarded0.778σ6, discarded0.758σ7, discarded0.658σ8, discarded0.571σ92.63·10⁻¹⁶σ102.3·10⁻¹⁶the operation that is not closedrank of each term4rank of the sum8truncated back to4cost of the truncation0.62the best there is0.62two planesmake a four-space
Fig. 24 The hierarchy field’s version of the same closure failure, from the essay that measures it.
How far the two residuals drift apart, against how large the iterates got on the waySix runs, differing only in the smallest eigenvalue of the matrix, which sets how far the iterates travel: from 9309 to 9.31·10¹³, ten decades of it. The upper line is the standard bound, u·‖A‖·max‖x‖, which every run stays inside. The lower series is what the drift actually does, and its fitted slope is 0.507 where the bound's is 1 — so the roundings are a random walk rather than a sum, the bound is correct, and it is loose by a square root that reaches 7.5·10⁴ across the sweep.10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.51a bound of one, a walk of a halffitted slope0.51the bound's slope1smallest gap measured5.3·10⁻¹⁵largest gap measured7·10⁻¹⁰share of the bound, near end7·10⁻⁴share of the bound, far end9.3·10⁻⁹the bound is rightand loose by a square root
Fig. 25 And the mechanism behind a residual that is not what a method thinks it is.
What residual replacement costs and what it buys, at three periodsEach row is the same conjugate gradient run with the recomputed residual assigned back into the recurrence every k steps. Without it, the reported residual reaches 4.01·10⁻²¹ and the answer's stalls at 6.24·10⁻¹¹. Replacing every 5 steps costs 30 extra matrix–vector products on top of 150 — 20 per cent — and brings the answer's residual to 7.17·10⁻¹⁷, with the two residuals then agreeing to a factor of 2.2.never replace6.24·10⁻¹¹replace every 252.02·10⁻¹⁶replace every 101.3·10⁻¹⁶replace every 57.17·10⁻¹⁷reported 4.01·10⁻²¹ · 0 extra productsreported 1.94·10⁻¹⁶ · 6 extra productsreported 8.22·10⁻¹⁷ · 15 extra productsreported 3.3·10⁻¹⁷ · 30 extra productsthe residual of the answer the run returnsone line, at three pricesnever: the answer's residual6.2·10⁻¹¹never: what it reported4·10⁻²¹every 5: the answer's7.2·10⁻¹⁷every 5: what it reported3.3·10⁻¹⁷extra products for that30ask the matrix againand the recurrence forgets what it did
Fig. 26 The repair that field found, which is available here and costs a dense operator application.
The residual of a formatted Cholesky, against the accuracy its blocks were compressed atA leaf of 32 means 10 truncations inside the factorisation, and the ranks of its blocks run to 4, 7, 9, 12, 14 across the sweep. The residual tracks the tolerance at a slope of 1.028 over ten decades: the knob does what a knob should, and nothing in the approximate arithmetic bends it. The ratio to the representation's own error stays at 0.94, 0.89, 0.92, 0.92, 0.93 — one excursion above one, at 0.01, where the tolerance sits between two singular values of one block and the rank it takes is a rounding of a decision rather than a measurement.10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²ε the blocks were compressed at‖A − LLᵀ‖ ⁄ ‖A‖the factorisationthe knob, on a factorisationtruncations10residual at 10⁻⁴1.1·10⁻⁵residual at 10⁻¹⁰9.6·10⁻¹²slope1worst ratio to the representation0.94ten decadesand a slope of one
Fig. 27 And what a tolerance buys inside a compressed factorisation, for reading against the rank budget here.

Where this goes next

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 The format that does not notice the dimension A Tucker core is r^d numbers, so the format that repaired the definition still cannot go past five indices. Cutting between the indices rather than across them gives d − 1 ranks instead of d, storage linear in the number of indices, and a family whose ranks are two everywhere by an addition formula. 02040608010010⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations2.1·10⁻¹⁰residual, 98 truncations1.1·10⁻⁹representation, deepest1.4·10⁻⁹residual ⁄ representation0.81levels, deepest5a hundred approximate stepsand an exact-looking factorisation The rounding that was not the problem A rank-k block plus a rank-k block is a rank-2k block, exactly, so every arithmetic in this format truncates after every addition. A Cholesky performed inside it does ninety-eight of those and its residual is 1.14·10⁻⁹ against a representation error of 1.40·10⁻⁹ — the roundings cost nothing measurable. 01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it The residual the method reports Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so.

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.

Conjugate gradientsCurse of dimensionalityKronecker sumKrylov subspaceLow-rank approximationModel problemRecompressionResidualTensor trainTruncation