When the problem arrives again

What survives one step of the barrier

An interior-point method solves the same system dozens of times with the same pattern and different numbers, and exactly p entries change between one step and the next. The pattern is reusable for ever. The factorisation is reusable for none of them, and the threshold that says so is a reduction factor of about a per cent against schedules that use ten.

Worth reading first: A factorisation kept past its date · The problem that arrives again · The zero that is not a missing entry.

This anchor’s first essay asked what a factorisation is worth when the matrix has moved a little, and answered it with a drift: a matrix perturbed by 10⁻³ relatively has a factorisation good for many steps of an outer iteration, and where the perturbation lands matters more than how large it is. The essay after it turned that into a rebuild policy with a crossover in it.

An interior-point method is the cleanest sequence of solves this collection has, and it gives the anchor’s question an answer with a threshold rather than a curve.

What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 6 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 1.02 and 2.88·10¹¹. The reason is in the entry count: the 6 entries that moved are the 6 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 1 What a factorisation survives, against the factor by which the barrier parameter falls. The number beside each mark is how many steps it lasted.

The sequence, and what changes in it

An interior-point method for a quadratic programme solves one saddle-point system per iteration:

[ H Cᵀ ] [ C −D⁻¹ ]

with D = Z/S the diagonal that the constraint field’s fifth essay is about. From one iteration to the next, H does not change, C does not change, and the sparsity pattern does not change. Exactly p entries move — the diagonal of the (2, 2) block, one per constraint — and no off-diagonal entry moves at all. Measured, at every step: 6 of 196 entries, 0 of them off the diagonal.

That is as close to “the same matrix again” as a sequence gets, and it is the case the anchor’s machinery was built for. The obvious expectation is that a factorisation of one step’s matrix is an excellent preconditioner for the next.

Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 3 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 29.04 to 3.044·10¹⁵ for the augmented form and 155.8 to 2.403·10¹⁶ for the condensed one, within a factor of 7.9 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.3098 and the augmented form's stays at 1.035·10⁻¹⁵ — fifteen correct digits at a condition number of 3.04·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable
Fig. 2 The matrix in question, from the constraint field, with the two ways of writing it.

What survives for nothing: the pattern

The sparsity pattern is identical at every step, so the symbolic factorisation is computed once for the whole optimisation and is free thereafter — provided the numeric phase never moves a pivot, which for the regularised form it never does. The sparsity field’s essay on that has the argument and the fill counts.

That is a genuine and large saving, and it is not what this page is about, because it is not a numerical reuse. Nothing about the values is being carried; what is carried is a shape.

The second thing that survives for nothing is the active set, once it has emerged. By μ = 10⁻⁸ the diagonal’s two groups are sixteen orders apart and no later step reclassifies anything. A code can use that: the split it reads off at one step is the split at every subsequent step, and that is what warm-starting an interior-point method means.

Fill under three orderings of a sparse saddle-point matrix, predicted symbolically and then countedA stiffness block of 24 unknowns with 4 constraints each touching 3 of them. Each panel is the pattern of PKPᵀ under one ordering, with the entries the factorisation adds shown against the entries the matrix arrived with. The counts are natural 113, minimum-degree 63, reverse-cuthill-mckee 63: computed from the graph before any number was touched, and matched exactly by the factorisation of the regularised matrix, because quasi-definiteness means nothing during the numeric phase is allowed to move a pivot. That is the guarantee sparse LU cannot make — there the pivot order depends on the values, so the symbolic phase can only bound the fill and the allocation has to be able to grow.natural113 predicted · 113 countedminimum-degree63 predicted · 63 countedreverse Cuthill–McKee63 predicted · 63 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural113minimum-degree63reverse-cuthill-mckee63predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 3 The pattern that is reused, and what it costs to allocate for.

What survives numerically: nothing

Carry the factorisation of K(μ) to K(σμ) and clean up with iterative refinement — the residual formed against the new matrix, the correction solved with the old factors. Twelve steps, at several reduction factors:

σ steps the factorisation survives 0.999 6 of 6 measured 0.99 6 of 6 0.95 2 0.9 1 0.5 0 0.1 0

Reuse holds while σ is within about a per cent of one, and fails outright by σ = 0.9. Real interior-point schedules reduce μ by a factor of two to ten, and at σ = 0.5 the first reused step already leaves a relative residual of 1.0 — no digits at all — while at σ = 0.1 it leaves 3.8·10⁹.

So the answer to the anchor’s question, on this family, is: every step. Not “rebuild when the drift exceeds a threshold” and not “rebuild every k steps”; rebuild every time, for a reason that is structural rather than a matter of tuning.

What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 4 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.0527 and 1.49·10¹⁰. The reason is in the entry count: the 4 entries that moved are the 4 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved4off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 4 At four constraints, where the threshold is in the same place — it is not a property of how many entries moved.

Why the intuition fails, in one number

Six entries out of a hundred and ninety-six sounds like a small update, and it is not, because those six are the six that dominate the matrix.

The relative change ‖K(σμ) − K(μ)‖ / ‖K(μ)‖ at σ = 0.1 is 9 — the matrix has changed by nine times its own norm in one step. The entries that moved are the ones of size 1/μ, and they moved by a factor of 1/σ, so the change is (1/σ − 1) times the largest thing in the matrix.

Put that on the same axis as the drift the anchor’s first essay measured, where a factorisation was still useful at a relative drift of 10⁻³ and marginal at 10⁻¹. A barrier step is a drift of nine. The reason the reuse fails is not subtle and not about where the perturbation landed: the perturbation is larger than the matrix.

The refusal this page publishes is the reading that counts entries instead of measuring them: a matrix that changes in only p of its entries between two solves is a small update, so a factorisation of the first is a good preconditioner for the second. It is fed the run at σ = 0.5 and required to fail.

Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 38 and 9 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁴.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10112233445566relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²38large end at 10⁻²9κ(M⁻¹A), small end36κ(M⁻¹A), large end1.3‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 5 The anchor’s original measurement, on the axis this one extends: where a drift lands, and what it is worth.

What the entry count is measuring, and what it is not

It is worth separating two readings of “six entries changed”, because the second is the one that makes the intuition wrong.

As a statement about structure it is exactly right and it is what the previous section uses: the pattern is unchanged, the elimination tree is unchanged, the allocation is unchanged, and everything a symbolic phase computed is still valid. Structural sparsity of an update is what lets a code avoid re-analysing.

As a statement about magnitude it says nothing at all, and the collection has met that distinction before from the other side. A rank-one update changes every entry of a matrix and can be applied to a factorisation in O(n²) because it is small in rank; a diagonal update changes p entries and cannot be applied at all when those p entries are the matrix’s norm. Sparsity of an update and smallness of an update are unrelated properties, and only the second governs whether a stale factorisation is any use.

The clean way to say it is that the relevant quantity is ‖ΔK‖·‖K⁻¹‖, and neither factor is a count. Here it is 9 times something of order one, and the count of six was never going to reveal that.

Sherman–Morrison against a direct solve, on 20×20 systems whose answer is the identityA is QΛQᵀ with Λ = (1, …, 1, ε) and the rank-one update takes ε back to 1, so A + uvᵀ is the identity and κ of the problem being solved is 1.000000 at every point on the axis. A direct solve returns 1.1·10⁻¹⁶. The update formula, which is exact algebra, returns 2.5·10⁻⁴ — a slope of 1.00 against κ(A), the matrix that was replaced. Its arithmetic is three solves against that matrix and there is nowhere else the error could have come from.10¹10³10⁵10⁷10⁹10¹¹10¹³10¹⁵10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A), the matrix that was updatedrelative forward errorSherman–Morrisondirect solve of A + uvᵀone answer, two routesκ of the answer's matrix1κ of the matrix replaced10·10¹³update formula's error2.5·10⁻⁴direct solve's error1.1·10⁻¹⁶the question's condition number is 1at every point on this axis
Fig. 6 The opposite case, from the least-squares field: every entry changed and the update cheap anyway.

The threshold is a property of the schedule, not of the problem

σ is chosen by the optimisation, for reasons about the outer convergence: an aggressive σ reaches the solution in fewer iterations and risks leaving the region where the Newton step is good, and every implementation has a rule for adapting it. Nothing in that rule is about linear algebra.

The measurement says the linear algebra’s reuse horizon is a function of exactly that number, and that the number the optimisation wants is two orders away from the number reuse would need. The counts do not move when the constraint count changes — the threshold is at σ ≈ 0.99 at four constraints and at ten — so it is not a property of the problem’s size either.

That makes it a clean example of something this field keeps finding: a policy question that looks like it should be tuned turns out to be settled, and settled by a quantity chosen elsewhere in the program for unrelated reasons. The rebuild policy essay found the optimum walking with a setup cost; here there is no optimum to walk, because the whole reusable region is on the wrong side of where the method operates.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 16.42 MFlop; refactorising every 5 costs 9.27. Past a period of 20 the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 4 factorisations, 9.22 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencedoes not convergethe contraction ruledrift 0.01 a memberevery member1.6·10⁷every 5 members9.3·10⁶contraction rule9.2·10⁶its factorisations4cliff at a period of20a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 7 The policy question in the case where it does have an answer, from this anchor’s own essay.

What refinement is doing, and why more of it does not help

The refinement iteration multiplies the error by roughly ‖I − K̂⁻¹K‖ each step, where K̂ is the stale matrix. That is a contraction only while the quantity is under one, and here it is approximately 1/σ − 1: about 0.01 at σ = 0.99, and 9 at σ = 0.1.

So the failure is a divergence rather than a slow convergence, and adding steps makes it worse rather than better. Twelve steps at σ = 0.1 leave a residual of 3.8·10⁹, and a hundred would leave more.

This is the same mechanism as the regularisation cliff in the constraint field, where refinement against the unregularised matrix works below δ = σ_min(K) and diverges above it. In both cases the useful reading is that iterative refinement is a fixed-point iteration and its convergence has a condition, and the condition is checkable before the iteration is run.

A code that reused a factorisation and refined until the residual stopped falling would get a diverging sequence and would have to detect it. A code that computed 1/σ − 1 first would not attempt it.

Iterative refinement on a factorisation of a matrix that is not the one being solvedThe last member of the sequence, solved with the first member's pivot order. Every other repair in this collection works because the outer loop recomputes the residual from the matrix, and this one does that too — and it does not work. The perturbed pivots are not a rounding, they are a change to the matrix, and the residual that would reveal it is computed with the same badly scaled rows that caused it. Six steps move the backward error from 4.82·10⁻⁹ to 5.46·10⁻¹⁰, a factor of 8.8, and stop. The equilibrated run has nothing to repair: it starts at 6.55·10⁻¹⁷.0123456710⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷refinement stepsbackward error of the solvethe kept order, in the units as giventhe same order, rows equilibratedthe repair that does not repairbefore refinement4.8·10⁻⁹after six steps5.5·10⁻¹⁰the factor gained8.8equilibrated, before6.5·10⁻¹⁷pivots replaced2the residual is recomputedand it cannot see what was done to the matrix
Fig. 8 Refinement failing for a different reason, from the field that owns it.

The one reuse that would work, and what it would cost

There is a way to make the factorisation reusable, and it is worth naming because it is not available and the reason is instructive.

The change between steps is D⁻¹(σμ) − D⁻¹(μ), a diagonal matrix supported on p entries. A factorisation can be updated for a diagonal change of rank p at a cost of O(p) rank-one updates, each O((n + p)²) for a dense factor or O(fill) for a sparse one — the machinery the least-squares field builds for exactly this shape. For p small against n that is cheaper than refactorising.

The obstacle is stability rather than cost. A rank-one update to a symmetric indefinite factorisation with a large multiplier is not backward stable — the update can be arbitrarily ill conditioned when the modification changes the inertia or nearly does, and here the entries being modified are the ones with the largest magnitude. The downdate essay in the least-squares field measures the same failure on a Cholesky, where removing an observation is a rank-one downdate that can lose every digit.

So the update exists, is cheap, and cannot be trusted on this family, which leaves refactorising. That is a common enough ending for an update argument that it is worth having the reason on the page rather than leaving the technique looking overlooked.

Downdating a Cholesky factor against refactorising it, as the removed row's leverage approaches oneA rank-one term is removed from a 6×6 Gram matrix by hyperbolic rotations and, separately, by factorising the downdated matrix from scratch. The downdate's residual rises from 2.1·10⁻¹⁶ at h = 0.3 to 3.5·10⁻¹⁰ at h = 1 − 10⁻⁷, a slope of 1.04 against 1/(1 − h). The refactorisation is flat at 8.1·10⁻¹⁷, and the matrix both are producing has a condition number of 4.3 at every point on the axis — so the difficulty belongs to the route and not to the answer.110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix4.3κ of the matrix downdated9.3·10⁶rotation's amplification344downdate residual3.5·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for
Fig. 9 The instability of the update that would otherwise work, from the field that measured it.

What could be reused instead, and is

The negative result is not the end of the matter, because there are three things a code does carry across a barrier step and all of them are cheaper than a factorisation.

The ordering and the symbolic factors, as above, which is most of the analysis cost.

The scaling. Equilibrating the matrix — which this collection has measured as making a reuse free elsewhere — depends on H and C, which do not change. So the row and column scalings are computed once.

The starting point. The previous step’s solution is a warm start for the next, and the whole design of a path-following method is that consecutive solutions are close. That is a reuse of the answer rather than of the machinery, and it is the one that pays: the outer iteration count is what it is because each Newton step starts near its target.

So the sequence is heavily reused and none of the reuse is a factorisation. Reading the anchor’s question as “how often to rebuild” gets a trivial answer here; reading it as “what is carried” gets three.

One sketch reused across 20 independent matrices, against a fresh sketch for eachThe matrices are drawn from the same construction with different seeds, so none of them can depend on the sketch. This is the case the randomised bounds were written for, and the two rows of points are the answer: mean error 0.36184 with one sketch kept and 0.35789 with 20 fresh ones, a difference of 0.00396 against a standard error of 0.00508. The rule is therefore not that a sketch must be redrawn — it is that the randomness is spent the moment the input starts depending on it.31323334353637383940410123relative approximation error ‖A − QQᵀA‖ ⁄ ‖A‖, per centone sketch, kepta fresh sketch eachthe vertical marks are the two meansreuse, where nothing adaptskept sketch, mean error0.36fresh sketches, mean error0.36difference0.004standard error of either0.0051spread, kept0.019the guarantee is about a drawthe input has not seen
Fig. 10 The anchor’s usual reading of what a carried factorisation is worth.

What an iterative method would do instead

Refinement is a stationary iteration and diverges when its contraction factor exceeds one. A Krylov method has no such requirement: it does not need the preconditioner to be close, only for the preconditioned spectrum to be clustered, and it will make progress from a preconditioner that refinement cannot use at all.

So a code that wants to reuse a stale factorisation across a barrier step has a route: use it as a preconditioner inside GMRES or MINRES rather than as a fixed-point map. The stale factors then buy a reduction in iteration count instead of an exact solve, and the question becomes how many iterations, which is a smoother trade than the cliff this page measures.

What the measurement here says about that route is not encouraging but it is quantitative. The preconditioned operator is K̂⁻¹K = I + K̂⁻¹ΔK, and ‖K̂⁻¹ΔK‖ is the same 9 that made refinement diverge, so the preconditioned spectrum is spread over a region of radius nine about one — not a cluster, and not obviously better than no preconditioner at all. At σ = 0.9 the radius is 0.11 and the clustering is excellent, which is the same threshold arriving in a different currency.

That is worth stating as a boundary on the negative result: reuse fails here for refinement and for a Krylov method, and it fails at approximately the same σ, because both are governed by the size of the same product.

LSQR steps to a normal residual of 10⁻¹², with and without a sketched preconditionerTwo curves against κ(A). The unpreconditioned solve takes 15, 22, 32, 59, 150 steps as the conditioning worsens, and at the worst one it does not reach the tolerance inside 150. The preconditioned solve takes 10, 10, 10, 10, 10 — flat across eight decades. Both return the same vector, to the accuracy the problem allows: the preconditioned answer and the direct one differ by about κ(A)·u at every κ, 5.7·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 4.5 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned10steps, without150answer gap, preconditioned5.7·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve4.5one sketch, two things to do with itand only one of them moves the answer
Fig. 11 Iteration counts under a preconditioner, in the field that measures what clustering is worth.

A control, and a boundary

The measurement needs a case where reuse does work, or the conclusion could be about the refinement machinery rather than about the barrier. σ = 0.999 is that case: six steps out of six at the rounding level, on the identical code path.

And the boundary between working and not is sharp enough to read off the figure — 0.99 works, 0.95 gives two steps, 0.9 gives one, 0.5 gives none — which is what makes the negative result a measurement rather than an absence.

The library’s other refusal covers the second half of the argument. It is fed the claim that a saddle-point matrix’s off-diagonal entries move between steps, and required to reject it: zero of them do, at every step, which is what makes the count of six meaningful and the drift of nine surprising.

What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 10 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.983 and 2.78·10¹¹. The reason is in the entry count: the 10 entries that moved are the 10 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved10off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 12 At ten constraints, where ten entries move per step and the threshold has not moved.

What this says about the anchor

Three essays into reuse, the shape of the answer has settled into a pattern worth stating.

Whether a factorisation can be carried is not a question about how many entries changed. It is a question about how much of the matrix’s norm changed, and about whether the change landed where the inverse is sensitive. The first essay found the second half of that; this one supplies a case where the first half is decisive and the entry count is actively misleading.

The interior-point sequence is the extreme instance: the sparsest possible update, in the densest possible sense. Six numbers out of two hundred, and a relative change of nine.

Chord steps per member of a drifting sequence, with the factorisation rebuilt every 4Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 3 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 12 at its worst, for 5 factorisations against 20.03691215182103691215member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 4factorisations5chord steps160worst member's steps12multiplications1.2·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 13 The anchor’s picture of a factorisation ageing, on a family where it ages gradually.

At other settings

What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 5 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.522 and 1.48·10¹¹. The reason is in the entry count: the 5 entries that moved are the 5 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved5off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 14 Five constraints, for reading against the hero’s six.
What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 8 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.915 and 2.59·10¹¹. The reason is in the entry count: the 8 entries that moved are the 8 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved8off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 15 And eight.
Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 1 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 28.23 to 1.439·10¹⁶ for the augmented form and 162.7 to 4.182·10¹⁷ for the condensed one, within a factor of 29.1 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 1 and the augmented form's stays at 1.352·10⁻¹⁵ — fifteen correct digits at a condition number of 1.44·10¹⁶.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴1.4·10¹⁶its relative error1.4·10⁻¹⁵κ₂ condensed4.2·10¹⁷its relative error1the same step, written two waysand only one of them is solvable
Fig. 16 The system, at one active constraint.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 3.044·10¹³ and 13.25 at μ = 10⁻¹², a ratio of 2.296·10¹², and the measured error is 9.434·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 13.25 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 1.324·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 17 And the reason its ill-conditioning is not what makes the reuse fail.
Relative error of a sum against the number of terms, at 24 significand bitsFive curves of relative error against the number of terms, both axes logarithmic. Three rise with unit slope; two rise with about half that.10²10³10⁴10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴terms summedrelative error+∞, −∞ 1.01zero 1.00stochastic 0.50nearest 0.47√n against n, fittednearest, fitted exponent0.47stochastic, fitted exponent0.5toward +∞, fitted exponent1twelve seeds averaged at each sizethe slope is the bias, not the precision
Fig. 18 The anchor’s drift measurement at a larger perturbation.
An incomplete Cholesky kept while the operator turns anisotropic, against one rebuilt at every memberThe five-point operator's two directions stop being equally weighted, ε running from 1 down to 0.0331, with the sparsity pattern unchanged throughout — so the factorisation from the first member stays applicable for the whole run, which is the situation in which it gets kept. The kept one goes from 18 iterations to 52. The rebuilt one goes from 18 to 13, because an anisotropic operator is an easier problem for a factorisation that knows about the anisotropy. The two start at the same point by construction and never meet again.0246810121416182001020304050member of the sequencepreconditioned conjugate gradient iterationskept from the first memberrebuilt every memberone operator, two policieskept, first member18kept, last member52rebuilt, first member18rebuilt, last member13ε at the last member0.033the problem got easierand the kept preconditioner got worse at it
Fig. 19 An incomplete factorisation ageing across a sequence, from this field’s own essay.
Iterative refinement of the inverse-and-multiply solve at κ = 10^8The backward error starts at 3.3·10⁻¹¹ and falls by about κu a step — 1.8·10⁻¹⁷, 2.7·10⁻¹⁷, 2.4·10⁻¹⁷, 2.1·10⁻¹⁷, 2.1·10⁻¹⁷ — reaching the LU route's 3.1·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.6·10⁻⁸ to 1.2·10⁻⁹, against the LU route's 9.6·10⁻¹⁰. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 3.1·10⁻¹⁷LU route: forward 9.6·10⁻¹⁰forward errorbackward errorwhat a correction buysη before refinement3.3·10⁻¹¹η after four steps2.1·10⁻¹⁷forward, unchanged1.2·10⁻⁹cost of a step, flops2048the residual is repairableand the accuracy floor is the problem's
Fig. 20 Refinement working as intended.
What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 21 The cost model an outer loop is priced in.
500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — [[H + δI, Aᵀ], [A, −δI]] with δ = 10⁻⁶ — has an LDLᵀ factorisation with diagonal D for every symmetric permutation, so the ordering can be chosen for fill with no numerical veto at all. That is the count: 500 of 500 orderings factorise here, against 69.0 per cent for the same matrix with the zero block left where the problem put it. Each bar counts the orderings whose factorisation left a relative residual ‖PKPᵀ − LDLᵀ‖/‖K‖ in that decade, with the tallest holding 136 of the 500. They span 1.56·10⁻¹⁶ to 1.02·10⁻⁷ — 9 orders — and the growth factor across them runs from 1 to 6.43·10⁵ — which is 0.643/δ, an inverse law that holds at every δ this figure is drawn at. Existence is not stability, and the theorem says only the first.-18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁵growth × δ0.64the ordering is free to chooseand not free of consequence
Fig. 22 The property that makes the pattern reusable at all.
What the regularisation costs, and what 6 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement6 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶3.6·10⁻¹⁶refined at δ = 10⁻²0.62the perturbation is known exactlybecause the code chose it
Fig. 23 And the other place a refinement iteration has a convergence condition with a number in it.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·10¹⁰ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 21.13 at the start of the sweep and 21.13 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 5.314·10⁻⁶ against 5.788·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 24 The alternative to factorising the whole thing every step.
What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 6 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.424 and 1.2·10¹¹. The reason is in the entry count: the 6 entries that moved are the 6 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 25 A larger primal block, where the six moved entries are a smaller share still.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 2.822·10¹³ and 13.17 at μ = 10⁻¹², a ratio of 2.142·10¹², and the measured error is 1.205·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 13.17 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 1.081·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²2.8·10¹³componentwise, same matrix13their ratio2.1·10¹²measured error1.2·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 26 The matrix’s conditioning at eight constraints.
Fill under three orderings of a sparse saddle-point matrix, predicted symbolically and then countedA stiffness block of 24 unknowns with 4 constraints each touching 4 of them. Each panel is the pattern of PKPᵀ under one ordering, with the entries the factorisation adds shown against the entries the matrix arrived with. The counts are natural 115, minimum-degree 67, reverse-cuthill-mckee 67: computed from the graph before any number was touched, and matched exactly by the factorisation of the regularised matrix, because quasi-definiteness means nothing during the numeric phase is allowed to move a pivot. That is the guarantee sparse LU cannot make — there the pivot order depends on the values, so the symbolic phase can only bound the fill and the allocation has to be able to grow.natural115 predicted · 115 countedminimum-degree67 predicted · 67 countedreverse Cuthill–McKee67 predicted · 67 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural115minimum-degree67reverse-cuthill-mckee67predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 27 The structure that is carried, at a wider constraint stencil.
What the regularisation costs, and what 3 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement3 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶4.2·10⁻¹²refined at δ = 10⁻²0.76the perturbation is known exactlybecause the code chose it
Fig. 28 And the refinement machinery, working where its contraction factor is under one.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Active setBarrier parameterCondition numberInterior point methodIterative refinementPreconditioningSaddle-point systemsSequence of solvesSparse pivotingSymbolic factorisation