Balanced is not symmetric

An answer that is known measured a linear solve against the Hilbert matrix’s integer inverse and ended by ranking the evidence an error claim can rest on: an answer known exactly, then an answer constructed backwards, then a backward error that certifies the algorithm and says nothing about the answer. The top level was available because the Hilbert matrix’s inverse is a closed form. It had a flaw the essay’s own later measurements exposed. The matrix a program holds is not the Hilbert matrix — every entry $1/(i+j-1)$ is rounded on the way into a double — and the exact answer to a nearby problem found that storing it accounts for between a fifth and all of what the solve gets wrong. The closed form is the answer to a matrix nobody stored.

An eigenvalue problem can be put at the top level without that flaw. The Sylvester–Kac matrix of order $n$ has zeros on its diagonal, $1, 2, \dots, n-1$ above it and $n-1, \dots, 2, 1$ below it. Every entry is a small integer, so the double a program holds is the matrix, with nothing rounded. And its eigenvalues are the integers $-(n-1), -(n-3), \dots, n-3, n-1$ — Sylvester stated it in 1854, and Kac met the same matrix again in 1947 as the Ehrenfest urn model of diffusion, where its eigenvalues are the rates at which a disturbed urn relaxes. The answer is known, the question is stored exactly, and every digit a solver loses is the solver’s.

Two exact routes to the integers #

The integers are not taken on trust. The characteristic polynomial of a tridiagonal matrix with zero diagonal obeys a three-term recurrence — each principal minor is $-\lambda$ times the one before, minus the product of the new off-diagonal pair times the one before that — and in exact integer arithmetic that recurrence is evaluated at every integer from $-(n-1)$ to $n-1$. At orders 8, 17, 40 and 120 it is exactly zero at the $n$ integers of the stated spectrum and nonzero at every integer between them. No tolerance enters: the spectrum is those integers and nothing else.

Three routes to the spectrum, and how far each gets #

The first route hands the matrix to the algorithm every library runs on a nonsymmetric matrix — reduction to Hessenberg form, then the Francis double-shift QR iteration, applied implicitly as two shifts that are never formed describes, after the reduction that makes it affordable. At order 8 the largest error is $1.6\cdot 10^{-14}$, which is rounding. At 64 it is $1.0\cdot 10^{-6}$; at 96, $4.3\cdot 10^{-2}$; at 104, 0.38, by which point neighbouring integers two apart are barely told apart. At 112 the solver returns eighteen of its 112 eigenvalues as complex numbers, the largest imaginary part 3.7, and at 120 twenty-six, with imaginary parts up to 6.4. The exact spectrum is 120 real integers.

Read as the urn model, the damage has a meaning. The eigenvalues are relaxation rates, evenly spaced by two, and at order 104 an error of 0.38 is a fifth of the spacing: the computed rates are still in the right order, but the gaps between neighbouring modes — which set how long an intermediate state of the urn persists — are wrong by up to thirty-five per cent. At 112 some modes no longer relax at all in the computed answer; they oscillate, because a complex eigenvalue is an oscillation, and the process has none. A user who read the answer’s structure rather than its digits would conclude something false about the physics, and no warning would say so.

The third route symmetrises first. A tridiagonal matrix whose off-diagonal pairs have positive products is similar, by a diagonal matrix, to the symmetric tridiagonal with off-diagonals $\sqrt{k(n-k)}$, and a symmetric eigensolver on that matrix returns the integers to $2.0\cdot 10^{-14}$ at order 8, $7.5\cdot 10^{-13}$ at 64 and $3.4\cdot 10^{-12}$ at 120 — within about twice $n$ units of roundoff of the norm at every order. The square roots are rounded where the Kac matrix’s entries were not; it makes no difference, because a symmetric eigenvalue moves by no more than the perturbation, which symmetry is worth more than precision is about.

Between them sits the route the libraries take by default: balance first.

The prediction, tested against integers #

The explanation for the first route’s collapse is the standard one, and it can be checked here more directly than anywhere else. The condition number of an eigenvalue, $\kappa(\lambda) = \|x\|\,\|y\|/|y^{\mathsf T}x|$ for its right and left eigenvectors, is the factor by which a perturbation of the matrix moves it. A condition number for one eigenvalue introduced it on a four-by-four; on the Kac matrix it can be computed for every eigenvalue at every order, because the eigenvectors are the symmetrised matrix’s eigenvectors scaled by the symmetrising diagonal and its inverse. A backward stable eigensolver perturbs the matrix by a few units of roundoff times its norm, so the prediction for each eigenvalue’s error is $\kappa(\lambda)\,u\,\|A\|$.

At order 48 the prediction is an arch — condition numbers of 2 at the two ends of the spectrum and $6.6\cdot 10^5$ at $\pm 1$ — and the errors of the unbalanced route sit on it, the largest $3.0\cdot 10^{-9}$ at $-3$. Turn the dial and the arch rises by fifteen or so every eight orders; at 96 its top is $6.6\cdot 10^{12}$, and the largest error is $4.3\cdot 10^{-2}$ against a prediction of $7\cdot 10^{-2}$. The extreme eigenvalues stay at rounding level at every order, because they stay well conditioned. What collapses is the middle.

Pooled over every order up to the collapse and both nonsymmetric routes, 1,568 eigenvalues lie along the diagonal over sixteen decades of condition number. The ratio of error to prediction has a median of 0.51; a tenth of the eigenvalues come out more than twelve times better than predicted, which is the prediction being a worst case over perturbations the solver did not happen to make; and the largest ratio is 135, at order 72 — the error a backward perturbation of about twice $n$ units of roundoff would produce, which is within what a Hessenberg reduction and a few hundred Francis steps are entitled to.

That is a stronger statement than the usual one, and it is available only because the answer is integers. The usual check of perturbation theory compares a computed eigenvalue with the same eigenvalue computed more carefully, which is the circularity the Hilbert essay opened with. Here the theory is checked against an answer that cannot be wrong in the same direction, over sixteen orders of magnitude, and it holds to within a factor that stays near one.

Balancing buys a constant and keeps the rate #

Balancing is the answer a library gives to a matrix like this one. Before the Hessenberg reduction, Parlett and Reinsch’s algorithm scales row $i$ down and column $i$ up by a power of two, repeatedly, until every row’s off-diagonal sum matches its column’s to within a factor of two. It is a diagonal similarity, so it changes no eigenvalue, and because the scalings are powers of two it rounds nothing. Its purpose is to shrink the norm the solver’s backward error is measured against, and with it the error. On the Kac matrix, which is a diagonal similarity away from symmetric, it is natural to expect it to find that similarity.

It does not. Balancing divides the worst condition number by 7.5 at orders 24 and 32, by 44 at 40, and by between 49 and 62 from 48 to 120 — a real gain, worth one and a half to two digits, and the hero figure shows it as a parallel shift of the error curve. What it does not change is the slope. As given, the worst condition number grows by a factor of about fifteen every eight orders, from 2.5 at order 8 to $2.3\cdot 10^{16}$ at 120. Balanced, it grows at the same rate to within a quarter, from 1.4 to $3.7\cdot 10^{14}$. The balanced solver reaches the unbalanced one’s errors about twelve orders later — sixty is fifteen to the power one and a half — and then collapses the same way: at 112 its largest error is 0.24, and at 120 it returns fourteen complex eigenvalues.

The symmetrised matrix’s condition numbers are all exactly one, at every order. So the diagonal similarity that would remove the whole of the growth exists, is cheap, and is what balancing is looking for. Balancing does not find it.

Fifteen every eight orders is a binomial coefficient #

The rate itself has a closed form, and it is worth deriving because it says what any diagonal scaling is up against. Symmetrising the Kac matrix means choosing $d_{k+1}/d_k$ so that the pair at position $k$ becomes equal, and the pair there is $k+1$ above and $n-1-k$ below, so the diagonal needs $d_{k+1}^2/d_k^2 = (n-1-k)/(k+1)$. Multiplying those ratios from the top gives $d_k^2 = \binom{n-1}{k}$: the symmetrising diagonal is the square roots of a row of Pascal’s triangle. Its entries run from one at the ends to the square root of the central binomial coefficient in the middle, which is about $2^{(n-1)/2}$, and an eigenvector of the symmetric matrix that spreads across the whole diagonal is stretched by that much on one side and shrunk by that much on the other.

Measured, the worst condition number is that square root divided by $\sqrt{n}$, times a constant: the constant is 1.20 at order 8, 1.14 at 32, and 1.131 from 64 to 104, which is within three parts in a thousand of $2/\sqrt{\pi}$ — a limit the measurement suggests and does not prove. Either way the growth is $2^{n/2}$ less a power of $n$, which is sixteen every eight orders with the power’s share taken off, and that is the fifteen the figure shows.

The derivation also says what balancing would have to do to remove it. The ideal scaling spans a factor of about $2^{n/2}$ from the ends of the diagonal to its middle — $10^{17}$ at order 120 — and balancing’s scalings span a factor of 64, at order 64 and at order 120 alike, because matching each row’s sum to its column’s is satisfied long before the pairs match. A diagonal that is allowed to span $10^{17}$ is exactly what the symmetrising route uses, in square roots rather than powers of two, and it is the only one of the three routes whose condition numbers do not grow.

Sums matched, pairs not #

The reason is in what balancing equalises. A symmetric tridiagonal matrix has, for each $k$, the entry above the diagonal equal to the entry below it — $n-1$ conditions on the pairs. Balancing imposes a different set: row $i$’s off-diagonal sum equal to column $i$’s, which for a tridiagonal matrix is $|a_{i,i-1}| + |a_{i,i+1}| = |a_{i-1,i}| + |a_{i+1,i}|$. A matrix can meet every one of those and have no equal pair at all. As given, the order-64 Kac matrix’s pairs run from $1/63$ to $63$. Balanced, they run within a factor of fifteen of one, and the worst balanced pair grows with the order — three at order 8, nineteen at 80, twenty-nine at 120 — because matching sums leaves a drift from one end of the diagonal to the other that the rounding to powers of two cannot remove.

The condition numbers follow the pairs rather than the sums. An eigenvalue condition number of a diagonally symmetrisable matrix is set by how unequal the scaled eigenvector’s two sides are, and that is the product of the pair ratios along the diagonal: an imbalance of fifteen at each of many positions compounds to the $10^{14}$ the balanced matrix still has at order 120. Balancing pulls each ratio in by a factor of about four at the ends of the diagonal and leaves the compounding in place.

This is the sense in which balanced is not symmetric. The two are different fixed points of different conditions, and a matrix can be at the first and as far from the second as the Kac matrix is. The ratio of the worst balanced pair to one is visible in a single pass over the matrix, which makes it a cheap warning: a balanced tridiagonal whose off-diagonal pairs are still unequal by ten or more has condition numbers balancing did not reach.

Where the spectrum goes when it collapses #

At order 120 the computed spectrum is not noisy integers. The outer eigenvalues, which are well conditioned, are on their integers to many digits. The middle ones — twenty-six of them unbalanced — have left the real axis and sit on a closed curve around the origin, a flattened oval reaching 6.4 above and below it; balanced, fourteen have left, on a curve reaching only 1.0. A backward stable solver has returned the exact eigenvalues of a matrix within $u\|A\|$ of the Kac matrix, and the Kac matrix is so far from normal that matrices that close to it have eigenvalues anywhere inside that curve.

That is the pseudospectrum drawn by the solver’s own rounding, and it is the same non-normality that makes a spectral radius that grows first climb before it decays. The eigenvalues that are not there drew the same region from the resolvent norm; here it appears as the set of answers a correct algorithm is entitled to return, and the exact integers make it possible to say that every point off the axis is wrong rather than merely different.

What an exact answer adds #

Three things were measured here that no comparison between two floating-point computations could have established. That the symmetrised route is right, rather than merely more consistent: it agrees with integers to $3.4\cdot 10^{-12}$ at order 120. That first-order perturbation theory predicts each eigenvalue’s error to within a median factor of two, not merely bounds it. And that the complex pairs at order 112 are an error of the whole method on this matrix, since the exact spectrum was verified integral before any of it was computed.

The Hilbert matrix could not give the first of these cleanly, because its question was already rounded before the solver saw it. The Kac matrix can, and the difference is worth stating as a rule for choosing test problems: a ground truth is only as good as the match between the problem it answers and the problem the program holds. An integer matrix with a known spectrum meets that rule; a matrix of fractions with a known inverse meets it only after the fractions are cleared, as a count that comes out of a determinant does with a Laplacian.

What the Kac matrix does not show #

It is one family, and a peculiar one: tridiagonal, zero diagonal, diagonally symmetrisable, with a spectrum of evenly spaced integers. Balancing’s failure here is a failure on the family it is least suited to — a matrix whose imbalance is spread thinly along a long diagonal — and on a dense matrix with a few badly scaled rows it does what it was designed for. The eigensolver is a direct implementation of the Francis iteration rather than a library’s, with the deflation test and shift strategy of the algorithm the libraries actually run, though the Hessenberg reduction and double shift are the textbook ones and the first-order ratio of 0.51 says it is behaving as a backward stable solver should. And the symmetrising route requires knowing that a symmetrising diagonal exists, which for a general nonsymmetric tridiagonal matrix is the sign condition on the pair products and for a dense matrix is almost never true.

Still open: the alternation, a balancing that matches pairs, and the eigenvectors #

Neighbours that differ tenfold. At order 40 the unbalanced errors alternate: $1.4\cdot 10^{-10}$ at $-9$, $8.4\cdot 10^{-12}$ at $-7$, $1.6\cdot 10^{-10}$ at $-5$, while the condition numbers change by a fifth from one to the next. Every eigenvalue and its negative have identical errors. The symmetry is the matrix’s — the diagonal of alternating signs is an exact similarity from the Kac matrix to its own negative, which is why the spectrum is symmetric about zero — and the prediction with a sign is that the alternation is the Francis iteration deflating the spectrum in pairs from the ends, so that eigenvalues found in the same deflation share their backward error. Recording which step deflated each eigenvalue would settle it.

A balancing that matches pairs. For a tridiagonal matrix the symmetrising diagonal is a single pass, one square root a position. For a general matrix no diagonal symmetrises, but a scaling that minimises the largest pair ratio, rather than equalising sums, is a small linear programme in the logarithms of the scalings. Whether it beats Parlett–Reinsch on dense nonsymmetric matrices with graded structure, or only on the Kac family, is the measurement.

The eigenvectors. The Kac matrix’s eigenvectors are known too — Krawtchouk polynomials, with integer entries — so the angle between a computed and an exact eigenvector can be measured the same way. The prediction is that the eigenvectors of the middle eigenvalues are lost before the eigenvalues are, since an eigenvector’s sensitivity carries the gap to the neighbouring eigenvalue as well as the non-normality, and the gap here is only two.