The elimination the matrix does not need

Two numbers retire the Kronecker form of a Sylvester equation, and both of them are already in this collection. At n = 100 the coefficient matrix of AX + XB = C has a hundred million entries for a problem whose answer is ten thousand numbers, and eliminating it costs (2/3)n⁶ = 6.67·10¹¹ operations against the 6·10⁷ that Bartels and Stewart’s algorithm needs. The essay that wrote the equation down states both as properties of the object and moves on to the conditioning, which is the more interesting question and also the reason nobody went back to look.

Neither number is a property of the object. Both are properties of a procedure applied to it, and the procedure is dense Gaussian elimination on a matrix whose entries nobody inspected. The entry count prices storing n⁴ numbers, and the matrix does not have n⁴ numbers in it. The operation count prices eliminating them, and after the reduction that both routes perform there is nothing left to eliminate.

What follows is that comparison redone with the structure counted rather than assumed. With A and B upper triangular — which is what a Schur reduction leaves, and the reduction is the first step of the algorithm that is supposed to be the alternative — the coefficient matrix, I ⊗ A + Bᵀ ⊗ I, holds exactly n³ nonzeros in its n⁴ entries, at n = 4, 6, 8, 10 and 12, and exactly zero of them above the block diagonal at every one of those sizes. Its diagonal blocks are themselves upper triangular. A matrix that is block lower triangular with triangular diagonal blocks is one permutation away from triangular, and a triangular system is solved by substitution rather than by elimination.

So the honest statement of the cost is not 6.67·10¹¹ against 6·10⁷. It is 10⁶ multiplications against 10⁶ multiplications, because the substitution that runs down the permuted Kronecker matrix and the substitution inside Bartels and Stewart’s algorithm are the same arithmetic in the same order. The two routes are not a method and a better method. They are one method drawn in the n⁴ picture and in the n² picture.

The matrix is sparse before anything is done to it #

Take A and B as they arrive, dense and with no structure at all, and form I ⊗ A + Bᵀ ⊗ I. The first term is n copies of A down a block diagonal, which is n·n² nonzeros. The second is a full n×n grid of scalar multiples of the identity, which is n²·n nonzeros. They overlap only on the diagonal of the whole matrix, where there are n² positions. So the count is 2n³ − n², and that is what the measurement returns: 112, 396, 960, 1,900 and 3,312 nonzeros at n = 4, 6, 8, 10 and 12, against n⁴ entries of 256, 1,296, 4,096, 10,000 and 20,736.

One number in 2.3 at n = 4 and one in 6.3 at n = 12, and the density falls as 1/(2n) forever after. At n = 100 the object with a hundred million entries holds 1.99·10⁶ numbers — sixteen megabytes rather than eight hundred, in any format that stores a sparse matrix rather than a rectangle. The n⁴ figure counts positions, not numbers, and the difference is a factor of n/2 before anything clever has been done at all.

That is not yet an argument, and it is worth being explicit about why, because the obvious next step is the one that fails. Sparsity is a property of a matrix and not of its factors. Elimination fills a sparse matrix in, and the amount it fills in is what decides whether the sparsity was worth anything. Run the symbolic elimination on this matrix in its natural order and count the structural nonzeros of the factors: 220 of 256, 1,146 of 1,296, 3,704 of 4,096, 9,190 of 10,000 and 19,284 of 20,736. 85.9%, 88.4%, 90.4%, 91.9% and 93.0% of n⁴, rising with the size.

So the (2/3)n⁶ estimate is not an exaggeration if the matrix is eliminated. The factors really are dense, the sparsity really is destroyed, and no ordering that reduces fill is going to recover a factor of n³ from a matrix that arrives 93% full at n = 12. The escape is not sparsity. It has to be something about the shape, and the shape is what the reduction supplies.

After the reduction there is nothing to eliminate #

Bartels and Stewart’s first step is two real Schur reductions, A = QRQᵀ and B = ZSZᵀ, which cost about 25n³ each and amplify nothing because they are orthogonal. The form a real matrix can reach is where that step is priced, and the case for taking it at all is that a similarity which changes no answer can change what the answer costs. The reduction is normally described as the step that makes the substitution possible. It is also, and this is the thing that was never counted, the step that makes the coefficient matrix triangular.

With R and S upper triangular, I ⊗ R contributes n copies of a triangle — n·n(n+1)/2 nonzeros — and so does Sᵀ ⊗ I, since Sᵀ is lower triangular with n(n+1)/2 nonzero entries and each of them is spread over an identity. Subtract the n² positions where the two diagonals coincide and the count is n²(n+1) − n² = n³, exactly. Measured on dense upper triangular A and B at n = 4, 6, 8, 10 and 12, the matrix has 64, 216, 512, 1,000 and 1,728 nonzeros, which are n³ with no lower-order term and no approximation. The density is exactly 1/n.

Three further readings come out of the same sweep, and each one is a fact the entry count cannot see. Zero nonzeros above the block diagonal, at every size — the first term is block diagonal, the second is block lower triangular, and no entry of the sum can be anywhere else. Every diagonal block upper triangular, being R + sⱼⱼI. And therefore: reverse the order of the unknowns within each block — send the entry (i, j) of X to position jn + (n − 1 − i) instead of jn + i — and the matrix is lower triangular, with zero entries above the diagonal, checked at all five sizes.

That permutation is not an ingenious device. It is the order the substitution already runs in: columns of X left to right, and each column from its last row to its first. The reordering that triangularises the Kronecker matrix is the loop order of the algorithm that never forms it, written down as a permutation instead of as two nested loops.

It is also load-bearing rather than cosmetic, which the same sweep says. In the natural order that stacking the columns of X produces, the matrix is not triangular: each diagonal block is upper triangular, so each contributes n(n − 1)/2 entries above the main diagonal, and the whole matrix carries n²(n − 1)/2 of them — 24, 90, 224, 450 and 792 at the five sizes, measured rather than counted by hand. A solver handed the matrix in that order finds a general sparse system and eliminates it. Handed the same matrix with the rows and columns of each block reversed, it finds a triangle. The structure was there both times; only one of the two orderings lets anything see it.

And both reductions are needed, which is worth checking rather than assuming, because A and B enter the sum in different ways and only one of them enters transposed. Reduce A alone and leave B dense: the count becomes n²(n + 1)/2 + n³ − n², which is 88, 306, 736, 1,450 and 2,520 nonzeros at the five sizes, above n³ at every one of them and spread over the full grid of blocks. The block diagonal is triangular and the block structure is not, so no permutation of this kind exists. Half the reduction buys a smaller count and none of the shape.

The two counts are one count #

A triangular system is solved by substitution, and a substitution costs one multiply-subtract per off-diagonal nonzero and one division per diagonal entry. The permuted matrix has n³ nonzeros, of which n² are on the diagonal, so the solve costs n³ − n² multiply-subtracts and n² divisions: n³ operations, exactly.

Bartels and Stewart’s substitution is the same total by a different argument. Entry (i, j) is obtained from a sum over n − 1 − i terms above it in its own column and j terms to its left in its own row, plus one division. Summed over the n² entries that is n²(n − 1) + n² = n³. Both derivations are worth having because they are counts of different things — nonzeros of a matrix nobody built, and arithmetic in a loop nobody wrote a matrix for — and they agree to the last operation.

Instrumented rather than derived, the two routes report 64 and 64 operations at n = 4, then 216 and 216, 512 and 512, 1,000 and 1,000, and 1,728 and 1,728. Not the same order of magnitude: the same integer, at every size. And the answers they produce agree to 4.4·10⁻¹⁷, 3.5·10⁻¹⁷, 1.6·10⁻¹⁶, 1.4·10⁻¹⁶ and 1.1·10⁻¹⁶ relative — which is rounding, and which is what agreement between two orderings of one arithmetic looks like.

The agreement extends to the arithmetic’s failure as well as to its success. The divisor at entry (i, j) is rᵢᵢ + sⱼⱼ, the sum of one eigenvalue of A and one of B, and it is the same divisor whether it is read off the diagonal of the permuted Kronecker matrix or computed in the loop. When it is zero the equation has no unique solution and both routes divide by zero at the same entry; when it is small both routes amplify by the same amount at the same entry. A triangular matrix’s condition number is not its diagonal, so that is not the whole conditioning story — sep is — but the point stands that there is no arithmetic in one route that is absent from the other.

The claim being made is therefore stronger than the Kronecker route can be made competitive. It is that after the reduction the Kronecker route is Bartels and Stewart’s substitution, with the indices flattened. There is no crossover, no size at which one becomes preferable, and no trade to report, because there is only one algorithm.

What the correction is worth, in the numbers that were quoted #

The dismissal was made at n = 100, so the correction belongs there.

Memory. The quoted figure is 10⁸ entries, 0.80 gigabytes of doubles. The reduced matrix has 10⁶ nonzeros, which is eight megabytes. A factor of a hundred, which is n, and it is the only factor of n that survives this whole argument.

Arithmetic. The quoted figure is (2/3)n⁶ = 6.67·10¹¹ operations. The substitution is n³ = 10⁶ operations. A factor of 666,667, which is (2/3)n³, and it is the price of eliminating a matrix that is already triangular.

And the sparse figure alone would not have been enough, which is the reason the shape had to be counted rather than the nonzeros. Stored sparsely without any reduction, the same matrix at n = 100 holds 1.99·10⁶ numbers — sixteen megabytes, a fiftieth of the quoted 0.80 gigabytes — and it is still the wrong object to build, because the elimination that follows fills it to more than nine tenths of n⁴ and the memory comes back. Sparsity moved the storage figure by a factor of n/2 and moved the arithmetic figure by nothing at all. Triangularity moves both, and it is the reduction rather than the sparsity that supplies it.

And the comparison the figures draw is unaffected, which is the part that has to be said plainly. Bartels and Stewart’s algorithm costs about 60n³ = 6·10⁷ at n = 100, the substitution is 10⁶ of that, and the rest is the two Schur reductions and four matrix products. The Kronecker route needs those reductions too — they are what made its matrix triangular. So the ratio of 11,111 between the curves is still the right comparison between forming and eliminating and reducing and substituting. What the measurement removes is the belief that those are two different algorithms rather than one algorithm and one unnecessary step.

Where it stops being true #

Three things, and the first is the one that matters in practice.

A complex pair breaks the triangle. The real Schur form of a matrix with complex eigenvalues is quasi-triangular: it carries 2×2 blocks on the diagonal, each with one entry below it, and the library’s own solver refuses such a case rather than pretending to handle it. Tensored with an identity, each of those subdiagonal entries becomes n entries of the coefficient matrix in the wrong place. Measured on random A and B at n = 6, 8 and 10, whose real Schur forms carried 4, 4 and 8 two-by-two blocks between them, the permuted matrix had 24, 32 and 80 entries above the diagonal — exactly n per block, and nothing else. The matrix is not triangular; it is triangular apart from a number of entries that grows like n times the number of complex pairs, which is why the fix is a 2×2 solve inside the loop rather than a different algorithm.

The factor of n is real and it is the reason nobody does this. n³ nonzeros against the n² + n numbers in R and S is a factor of about n, and against the n² entries of X it is a factor of exactly n. At n = 100 that is eight megabytes for a coefficient matrix beside eighty kilobytes for the two triangular factors it was built from. At n = 400 it is 512 megabytes beside 1.3. Storing the same information n times over is not fatal and it is not free, and it buys nothing, because the substitution reads R and S directly.

And the index arithmetic is not free either. A sparse triangular solve chases pointers; the substitution written in two loops reads contiguous rows. Both cost n³ arithmetic operations and the second is faster, which is the ordinary reason a structured algorithm beats a general one applied to a structured object — the same reason the two-index form of a grid is worth keeping rather than flattening.

What this does not license #

The correction says the object is cheaper than advertised. It does not say the object should be formed, and the distinction is the one this collection draws everywhere.

Forming it is still a factor of n in memory for no gain, which is the whole of the practical verdict. The site’s licence for building an object that ought not to exist is narrow and stated: form it when its entries or its singular values are the answer, never as a step towards one. sep is the smallest singular value of this matrix and there is no other way to say what it is, so measuring the conditioning of a Sylvester equation is exactly the case the licence covers. Solving the equation is exactly the case it does not.

And the refusal this page publishes still refuses. At n = 20 the elimination costs 4.27·10⁷ operations against 4.8·10⁵, and the assertion that the two are within a factor of ten is handed that comparison and fails on it. What the measurement changes is the reason it fails. It is not that the coefficient matrix is a bad object; it is that the flop count prices an elimination the matrix has no need of, and a claim that survives being priced correctly is a better claim than one that needed the overstatement.

The same shape appears wherever a small object is written as a large one. A rank-one update is a matrix that is one row and its n² entries are notation. A Kronecker sum of d factors is d decompositions rather than one, and the field next door solves the same object by diagonalising each factor instead of triangularising it. An intermediate that is never needed as a matrix is a product nobody had to form. In every case the entry count of the written-down object is a count of positions and the algorithm is about the numbers.

What is worth carrying #

A count of entries is not a count of numbers. The coefficient matrix of AX + XB = C has n⁴ entries and 2n³ − n² numbers before the reduction, and n³ after it. Quoting the first as the cost of the object is the error, and it is an easy one because n⁴ is what the notation says.

A count of operations is a count of operations in a stated algorithm. (2/3)n⁶ is the cost of dense elimination, and dense elimination is the wrong algorithm for a matrix that is a permutation away from triangular. The measured cost of the right one is n³, exactly, which is the substitution the standard method performs.

The reduction is the load-bearing step in both pictures. It is what makes the substitution possible and what makes the coefficient matrix triangular, and it is the same two orthogonal similarities either way. Anything said about the Kronecker form that does not say whether the reduction has happened is ambiguous between a matrix that is 93% full after elimination and one that needs no elimination at all.

And a crossover between two methods should be counted before it is quoted. The flop ratio here is n³/90, so it passes 1 at n = 4.48 and reads 0.71 at n = 4 and 2.4 at n = 6 — which is where the received advice about small n comes from, and it prices the wrong route on both sides of the crossing. Where a format starts paying is a question with an answer, and the answer is a measurement rather than an estimate. Here it turns out there is no crossing to find: the two curves being compared were never two algorithms, and the object drawn at the top of a two-dimensional model problem is the same Kronecker sum, arriving from a discretisation instead of from an equation.