An iteration that only multiplies

The polar factor of the previous essay was computed through a singular value decomposition, which is the most expensive factorisation in the subject. Nobody who wants to re-orthogonalise a rotation every millisecond is going to pay for one.

There are two iterations that get U without it, and the difference between them is the subject here. One converges faster and needs a factorisation every step. The other needs nothing but matrix products — no pivot, no comparison, no entry ever read — and pays for that in a step count and in a boundary.

Newton, and what it is Newton’s method for #

The orthogonal matrices are exactly the solutions of X⁻ᵀ = X. Applying Newton’s method to that equation gives

X ← (X + X⁻ᵀ)/2

started at A, and it converges to the orthogonal polar factor. It is a beautiful iteration — the scalar version is x ← (x + 1/x)/2, which is Heron’s method for √1 and converges to 1 from any positive start — and each step costs one LU factorisation and n triangular solve pairs, which is 2n³ operations.

Quadratic convergence, an inverse a step. And on a badly conditioned matrix it is much worse than that.

The first phase is linear, from a quadratic method #

At κ = 10⁶ the unscaled iteration’s error goes

1.8·10⁵ → 8.9·10⁴ → 4.5·10⁴ → 2.2·10⁴ → 1.1·10⁴ → 5.6·10³ …

falling by exactly 0.5000 a step. That is linear convergence, from a method advertised as quadratic, and the ratios are 0.49999886, 0.49999772, 0.49999544 — not approximately a half.

The scalar map explains it. Each singular value σ follows σ ← (σ + 1/σ)/2 independently. For σ ≫ 1 the 1/σ term is negligible and the map is σ ← σ/2: halving. Only once σ is within a factor of a few of 1 does the quadratic behaviour start.

So the number of steps is about log₂ σ_max plus a handful, and σ_max is whatever the matrix’s largest singular value happens to be. A method whose step count depends on a scalar multiple of its input is in the same position as the Taylor series two essays ago, and for a related reason.

The scaling, which is free #

Higham’s remedy multiplies X by ζ and X⁻ᵀ by 1/ζ before averaging, with

ζ  =  ( ‖X⁻¹‖₁ ‖X⁻¹‖_∞ / ‖X‖₁ ‖X‖_∞ )^(1/4)

which is an estimate of 1/√(σ_max σ_min): it moves the singular values so that the largest and smallest are equally far from one. The inverse is already being computed, so the two extra norms are O(n²) on a step that costs O(n³) — the scaling is free.

Six steps against more than twelve, for two norms. This is worth stating as a habit rather than as a fact about one iteration: where a method’s early phase is governed by a scale, a scaling computed from quantities already in hand is almost always available and almost always ignored.

The iteration that never inverts #

Newton–Schulz is

X ← X(3I − XᵀX)/2

Two matrix products a step. No factorisation, no pivot, no triangular solve, no comparison of any entry with any other. On a machine where a matrix multiply is a hundred times more efficient than a factorisation — which is to say, on every machine built in the last twenty years — that is a very different proposition from 2n³ of triangular solves.

Its whole theory is a scalar cubic. Each singular value follows

σ ← σ(3 − σ²)/2

whose fixed points are 0, 1 and −1, and whose basin of attraction around 1 is exactly the interval (0, √3).

The boundary is sharp to five digits #

Scale A so that its largest singular value is exactly α, and run the iteration to convergence:

and √3 = 1.7320508. That is a two-sided measurement of a boundary, which is the form this site reaches for whenever a slider crosses one, and it is as sharp a boundary as the collection has.

What it does past the boundary is worse than diverging #

Here is the finding, and it is not the one the theory suggests.

At α = 1.73206 the iteration converges. It converges to a matrix whose orthogonality residual ‖XᵀX − I‖ is 5.1·10⁻¹⁶ — orthogonal to the last bit — and whose relative error against the polar factor is 0.707.

One of the singular values was carried to −1 instead of +1. The result is an orthogonal matrix with a reflection in it, and nothing in the output says so. A caller checking that the answer is orthogonal gets a yes. A caller checking that the iteration converged gets a yes. The determinant is −1 instead of +1, which would catch it — and the determinant is the scalar this phase opened by arguing against, so it is not what anybody computes.

And past that, the limit is not a function of anything #

The cubic’s behaviour outside (0, √3) is not a simple divergence. Three starting values agreeing to five digits:

2.23600  →   1
2.23606  →  −1
2.23607  →  ∞

The second boundary is √5 = 2.2360680. Between √3 and √5 the map folds the interval onto a wider one and the dynamics are chaotic: which fixed point a starting value reaches is not a continuous function of it, and no amount of care in choosing α helps.

Only the third of those three outcomes is visible from the output. The first is right, the second is orthogonal and wrong, and the third overflows.

So the practical rule is not check whether it converged. It is: scale first, and scale conservatively. Dividing A by ‖A‖_F guarantees every singular value is at most 1, which is safely inside the basin, and costs one pass over the matrix.

Why the residual is the wrong stopping test here #

Every iteration in this essay has to decide when to stop, and the obvious test is unavailable for an interesting reason.

For a linear solve the residual b − Ax̂ is computable and cheap, and it decides. For the polar factor the natural analogue is ‖XᵀX − I‖ — how orthogonal is the current iterate — and it does not decide, because it is small for every orthogonal matrix and the iteration has more than one orthogonal fixed point.

The measurement above makes that concrete: past the basin, ‖XᵀX − I‖ reaches 5·10⁻¹⁶ and the answer is wrong. A stopping test built on it stops, and reports success, on a matrix with a reflection in it.

What does decide is the step size — ‖X_{k+1} − X_k‖ — because the map is a contraction near its fixed point and the increment goes to zero only at one. It costs nothing beyond the step already being taken. It is also, notably, the test that a Newton iteration for a scalar root would use, and the reason the matrix version reaches for the residual instead is that the residual exists, which is a poor reason.

A quantity that is small for every fixed point cannot distinguish between them. That sentence is worth carrying past this essay: several of this site’s methods have more than one fixed point and the convention is to test the one thing they have in common.

What the inverse-free method charges #

Step counts to thirteen digits, on matrices of prescribed conditioning:

The slope is 5.67 steps a decade of κ, and log 10/log 1.5 = 5.68. That is not a fitted constant with a plausible story attached; it is the closed form. A singular value below one is multiplied by about 3/2 per step while it is small — the map’s derivative at zero is 3/2 — so the number of steps to bring σ_min up to one is log κ / log 1.5.

So the trade is exact and stateable: the inverse-free method pays for what it does not do in a step count proportional to log κ, at two matrix products a step, against a method that pays a factorisation a step and needs six of them.

Which is better depends on the ratio of a matrix multiply to a factorisation on the machine in hand, and on κ. At κ = 10 and a multiply that is fifty times cheaper, Newton–Schulz wins comfortably. At κ = 10⁶ it needs about forty steps and the arithmetic is close.

Where the products go, which is not obvious #

Two matrix products a step sounds like a small thing until the sizes are put in.

Newton–Schulz on an n×n matrix is 2 × 2n³ = 4n³ operations a step, against a Newton step’s factorisation and n solve pairs, which is also about 2n³. So per step the inverse-free method is twice the arithmetic, not half — and it needs more steps.

The reason it is preferred anywhere is that the arithmetic is not the cost. A matrix multiply runs at close to a machine’s peak rate: it is O(n³) operations over O(n²) data, so every number loaded is used n times, and every processor built since the nineties has been designed around exactly that access pattern. A triangular solve is O(n²) operations over O(n²) data and is a sequence of dependent steps; it runs at a fraction of peak, and on a distributed machine it needs a global communication per column.

This is the cost field’s whole argument arriving in an orthogonalisation. The operation count says Newton wins by a factor of five or ten; the measured time can go either way, and on a GPU or a large distributed machine it usually goes the other way.

That is the honest reason to write down a method that does twice the arithmetic and takes four times the steps. It is not elegance and it is not accuracy; it is that a matrix multiply is the operation hardware is built for and a factorisation is not.

The same iteration under three names #

The Newton iteration for the polar factor, the Newton iteration for the matrix sign function, and the Newton iteration for a matrix square root are the same iteration wearing different clothes.

sign:    X ← (X + X⁻¹)/2
polar:   X ← (X + X⁻ᵀ)/2

and sign(A) applied to the block matrix [0 A; Aᵀ 0] returns [0 U; Uᵀ 0] with U the polar factor. One iteration, three literatures, three sets of convergence proofs written independently.

Newton–Schulz has the same triple life, and in each of them it is reached for under the same pressure: a setting where a factorisation is expensive or impossible. Lattice field theory applies the sign function to an operator of order 10⁷ that exists only as a subroutine. Machine learning orthogonalises weight matrices inside a training loop where a matrix multiply is what the hardware does. Neither could invert anything if it wanted to.

Higher orders, which buy steps with products #

Newton–Schulz is the first member of a family, and the family is worth knowing because it makes the trade explicit rather than fixed.

The scalar map σ ← σ(3 − σ²)/2 is the order-2 member of a sequence obtained by truncating the binomial series for (1 − x)^(−1/2). The next one is

X  ←  X (15I − 10 XᵀX + 3 (XᵀX)²) / 8

which converges cubically rather than quadratically, costs four matrix products a step instead of two, and has the same basin — the polynomials all fix ±1 and all have their real basin bounded by where the derivative condition fails. So a higher order buys a lower step count at a proportionally higher cost per step, and the product of the two is roughly constant near convergence.

Where it is not constant is at the start, on a badly conditioned matrix, and that is the half worth knowing. The convergence there is governed by how fast a small singular value is amplified: the order-2 map multiplies it by 3/2 a step and the order-3 map by 15/8, so the number of steps to bring σ_min up to one goes from log κ/log 1.5 to log κ/log 1.875, which is a factor of 1.55 fewer steps for a factor of 2 more products. The higher-order member is worse on the approach and better at the end, and the practical answer is to switch: run the cheap map until the iterate is inside a small neighbourhood of orthogonality, then finish with the expensive one.

That is a familiar shape — it is the same reason Newton’s method is preceded by a bisection, and the same reason the scaled Newton iteration exists at all. An iteration with a slow approach and a fast end is two algorithms with a hand-over, and the interesting parameter is where the hand-over goes rather than which iteration is better.

What is worth carrying #

A quadratically convergent method is linear until it is close, and close can be twenty steps away. The unscaled Newton iteration halves its error for the first phase, and the length of that phase is log₂ of the largest singular value.

A scaling built from quantities already computed is free and is worth six steps here. Two norms of a matrix whose inverse is being formed anyway.

And an inverse-free iteration has a basin, sharply. Below √3 it converges; above it, it converges to an orthogonal matrix that is not the answer and reports nothing; above √5 it overflows. Scale into the basin rather than detecting the failure, because two of the three failures are invisible.

The last essay of this phase takes the property that made Newton–Schulz attractive — that it never reads an entry — and asks what an operator with no entries at all costs: an operator with no entries.