A rotation that comes back mirrored

The nearest orthogonal matrix ends with a warning it does not measure. The polar factor minimises ‖A − Q‖ over every orthogonal matrix, and the orthogonal matrices come in two pieces: the rotations, with determinant +1, and the reflections, with determinant −1. An application fitting a rotation to two point sets — a molecule aligned to a reference, a spacecraft’s attitude from star sightings, a scan registered to a model — wants the first piece only, and the polar factor is free to return the second. The repair in every Kabsch implementation is to flip the sign of the singular vector belonging to the smallest singular value. “Correct on every input where the data is not nearly degenerate,” the warning runs, “and catastrophically wrong on the ones where it is.”

Three things in that sentence are checkable, and none was checked. How often does the polar factor come back mirrored, and as a function of what? How wrong is “catastrophically” — a large error, or a different kind of error? And nearly degenerate in which direction, since a point set can be degenerate by being flat or by being thin like a rod, and the fix might not treat both alike?

The experiment: a thin disc, turned and blurred #

Twenty points are drawn as (z₁, z₂, t·z₃) with each z standard normal: a disc of unit spread and thickness t. A rotation Q★ with determinant +1 is drawn at random, the points are rotated by it, and independent noise of standard deviation σ is added to every coordinate of the result. That gives two point sets, B and A = BQ★ + noise, and the problem is to recover Q★ from them.

The least-squares answer over orthogonal matrices is the polar factor of the 3×3 cross-covariance M = BᵀA. Write M = PΣQᵀ; the polar factor is U = PQᵀ. It is the nearest orthogonal matrix to M, and the nearest orthogonal matrix measures that it is nearer than QR’s answer and than two hundred random draws. The rotation with the same minimisation restricted to determinant +1 is R = P·diag(1, 1, d)·Qᵀ with d = det(PQᵀ). When d = +1 the two are the same matrix. When d = −1 the smallest singular direction is flipped.

The determinant here is the right test and a rare case where the scalar the number that decides nothing argues against is exactly what is wanted: det(PQᵀ) is ±1 by construction, so its sign is decided by a margin of one and there is no threshold to choose. It is computed by cofactors and checked against the LU product on every mirrored trial; the two never disagree.

How often: a rate set by noise over thickness #

At each ratio σ/t from 1 to 100, 400 alignments are run at each of three thicknesses, 10⁻², 10⁻³ and 10⁻⁴, each thickness on its own seed so that agreement between them is a measurement rather than the same random numbers scaled. The percentage of mirrored polar factors, averaged over the three:

The thickness does not matter; the ratio does. At σ/t = 3 the three thicknesses give 6.3, 9.3 and 7.5 per cent; at 10, 37.3, 31.5 and 33.8. Those spreads are what 400 trials a cell produce — the figure’s own check allows four standard deviations of the difference of two binomial proportions and none is close to it — and they come from point sets whose thicknesses differ by a factor of a hundred. A disc of thickness 10⁻⁴ measured with noise 10⁻³ is mirrored exactly as often as a disc of thickness 10⁻² measured with noise 10⁻¹.

And the rate rises towards a coin, not past it. At σ/t = 100 it is 46.6 per cent. Once the noise swamps the thin direction entirely, nothing in the data distinguishes a disc from its mirror image, and the polar factor picks one at random.

The mechanism is the smallest singular value of M. The thin direction contributes about m·t² to it from the points themselves, and the noise contributes a term of size about σ·t·√m with a random sign. When the second is comparable with the first, the smallest singular direction’s orientation is decided by the noise, and the determinant of PQᵀ with it. That heuristic predicts the ratio σ/(t·√m) as the governing quantity rather than σ/t, and it is a heuristic: the three-thickness agreement is checked, the √m is not.

The point count does move the curve, and in roughly the way the heuristic says. At σ/t = 10 the mirrored share falls from 43.0 per cent with five points to 38.9 with ten, 34.2 with twenty, 26.4 with forty and 18.4 with eighty. Eighty points at a ratio of ten mirror about as often as twenty points at five, and forty points at ten fall between twenty points at five and at ten. Five points do not follow the pattern closely — five points in three dimensions is a small sample, and the distribution of σ₃ is not yet near its large-m shape — which is one reason the dependence on m is reported here and not checked.

How large: exactly two, and exactly 4σ₃ #

“Catastrophically wrong” suggests a large error. The measurement says it is not an error of a size at all.

Of 300 alignments at σ/t = 10, 198 polar factors are rotations and 102 are reflections. For the 198, ‖U − R‖ is 0 — not small, zero, since U and R are the same product. For the 102, ‖U − R‖ is 2 to within 8.9·10⁻¹⁶. No trial falls anywhere between. The distance is a property of the flip rather than of the data: U − R = P·diag(0, 0, 2)·Qᵀ, whose Frobenius norm is 2 because P and Q are orthogonal.

So the polar factor of a noisy alignment is either the rotation exactly or the rotation composed with a half-turn reflection across a plane, and there is no intermediate state in which it is a little bit mirrored. A code that checks the distance from the previous frame’s U for smooth changes will see nothing, then a jump of 2, then nothing — the signature of a discrete event, not of accumulating error.

The residual says how strongly the data preferred the mirror, and the answer is an identity. Expanding ‖A − BX‖² for orthogonal X gives ‖A‖² + ‖B‖² − 2·trace(XᵀM), and with X written in the singular vectors of M the trace is σ₁ + σ₂ ± σ₃. So the rotation’s residual exceeds the reflection’s by exactly 4σ₃. In every one of the 102 mirrored trials that identity holds to a relative 8.3·10⁻¹³, which is two routes to one number — two subtracted residual norms on one side, one singular value on the other — agreeing to the precision the subtraction allows.

And σ₃ is small. Among the mirrored trials it runs from 3.6·10⁻⁵ to 7.1·10⁻³, and the median price of the flip is 1.1·10⁻² against a total squared residual of about 0.6. The data prefers the reflection by a couple of per cent of its own misfit. That is the quantitative content of “nearly degenerate”: not that the rotation is poorly determined, but that the handedness is decided by a residual difference the noise can manufacture.

What the fix buys: an error that ignores the thickness #

If the mirror is a discrete event, the interesting question is what the rotation’s error is once the event is removed — and whether a thinner disc, which carries less information about its third direction, gives a worse rotation.

Thinning the disc costs the rotation nothing. From a thickness of 0.03 down to 3·10⁻⁴ — two orders of magnitude, including the whole range where the polar factor is failing — the fixed rotation’s error stays between 0.476 and 0.513 of the noise. The one change in its column is between the ball and the disc: a ball has a third direction of spread that helps pin down the rotation, and losing it raises the error by about a third, from 0.375 to 0.49. After that, the thin direction carries no more information to lose.

That is worth a sentence of interpretation, because it is the reason the fix is a fix rather than a patch. The thin direction’s only job in the estimate is to decide the handedness, and the handedness is known before any data is taken: it is +1, because the thing being estimated is a rotation. Imposing det = +1 supplies exactly the information the thin disc stopped carrying. The two in-plane directions determine the rotation, the constraint determines the third axis’s sign, and the data’s thinness stops mattering.

The polar factor’s column is the same computation without the constraint, and it tracks the rotation to two figures while no trial is mirrored — 0.38, 0.47, 0.49 — then departs as soon as mirrors begin. Its error at thickness 3·10⁻⁴ is an average of 0.48σ over the unmirrored trials and 2 over the mirrored ones, which is why it reads as 90σ: a 2 divided by a noise of 0.01, diluted by the trials that were not flipped.

The threshold moves with the noise and not with the thickness alone, which is the first figure’s finding seen from the other axis. At σ = 10⁻³ no trial is mirrored down to a thickness of 10⁻³, and 11 per cent are at 3·10⁻⁴. The rotation’s error per unit of noise is the same column as before. The alignments are drawn with the same seeds at both noise levels, so the unchanged column says that the error is first-order in the noise over this range rather than confirming it on independent draws.

The fix has a condition number, and a rod pays it #

A correction that changes the answer by exactly 2 when it fires is discontinuous, and a reasonable worry is that the corrected rotation is badly conditioned near the switch. It is not — near the switch. It is badly conditioned somewhere else.

The shapes here are built rather than drawn: the cross-covariance is W·diag(1, s₂, ±s₂·q)·Wᵀ·Q★, so σ₃/σ₂ is exactly q and the mirrored shape has the same singular values as the unmirrored one. That matters because a random draw of point coordinates does not get close to σ₂ = σ₃ — sampling noise holds them apart by a factor of two — and the limit this measurement is about would never be reached.

On a disc, σ₃ ≈ 0 and the two formulas agree: flipping the smallest direction of a disc is a well-conditioned operation, 15.7 against 13.9. That is the case every figure above is about, and it is why the fixed rotation’s error did not grow as the disc thinned.

On a rod, σ₂ ≈ σ₃ and the flipped formula diverges: at σ₃/σ₂ = 0.98 the rotation moves 724 times the perturbation, against 7.4 without the flip. A rod is a point set spread along one axis with two nearly equal short axes, and the determinant constraint forces a choice of which short axis to reflect across — between two directions the data can barely tell apart. A small change in M swaps them, and the answer rotates by a large angle about the rod.

That instability is not an artefact of the fix. Rotation about the long axis of a rod is poorly determined by the data in any formulation, and the unflipped rotation’s sensitivity would grow as the short axes shrank too; what the flip adds is a factor of (σ₂ + σ₃)/(σ₂ − σ₃), which is harmless for a disc and unbounded for a rod whose short axes are equal. The measured values sit within a quarter of their formulas at every stop, and the figure requires agreement to within a factor of 1.6 in either direction — forty random perturbations find the worst direction only approximately.

So “nearly degenerate” has two meanings with opposite consequences. A flat point set makes the polar factor wrong often and the fixed rotation not at all. A rod makes the fixed rotation sensitive whenever the fix is needed. The warning that started this essay was about the first and is silent about the second.

What a caller can read off the data #

Every measurement above uses a truth that a real alignment does not have, so it is worth separating what a caller can compute from what only the construction knows. The determinant of PQᵀ is available and is decided with a margin of one, which makes it the kind of test an eigenvalue count that cannot be slightly wrong describes: the number itself cannot be slightly off. What it cannot say is whether the sign it reports belongs to the object or to the noise.

The quantity that says that is also available, because it is the price of the flip. The rotation’s residual exceeds the reflection’s by exactly 4σ₃, and σ₃ is the smallest singular value of the cross-covariance that the fix already computed. In the mirrored trials at σ/t = 10, 4σ₃ runs from about 10⁻⁴ to 3·10⁻², against a total squared residual of about 0.6. A caller who compares 4σ₃ with the spread the residual would show under a fresh draw of the noise is asking whether the data’s preference for one handedness is larger than the noise could manufacture — a comparison with a threshold in it, where the determinant had none.

That comparison is argued here and not measured. The measurement it would need is the distribution of 4σ₃ across trials that are and are not mirrored at a given ratio, and the fraction of each that a threshold on it would flag. The figures above give one end of that distribution — every mirrored trial at σ/t = 10 has σ₃ below 7.1·10⁻³ — and not the overlap with the unmirrored trials, which is the number that decides whether such a check is worth running.

The shape of the decision is familiar from elsewhere in the collection. Deciding that a zero has arrived is about detection thresholds where both kinds of error are failures, and a handedness check is one: flagging a correct rotation as doubtful costs a caller a recomputation, and passing a mirrored one costs them a molecule reflected through a plane. The fix itself removes the second failure whenever the physics says the answer is a rotation, which is why the essay’s recommendation is the fix and not the check. The check is for the cases above where the polar factor is still wanted, and a caller needs to know whether the handedness it returned is the data’s or the noise’s.

When the polar factor is still the answer #

None of this makes the polar factor wrong. It answers “the nearest orthogonal matrix”, and for several uses reflections are legitimate. Re-orthogonalising a drifted matrix that was orthogonal but not necessarily proper — a basis, an orthogonal transformation of coordinates that may include a handedness change — wants the unconstrained answer, and the flip would be the error. A test with no answer in it checks a routine by reversing its columns, and a column reversal of a 3×3 matrix changes the determinant’s sign: the polar factor passes that test because it permutes with the columns, and the fixed rotation, which must keep determinant +1, does not permute with an odd permutation at all. Which of the two a caller wants is a statement about the physics.

The iteration that avoids the SVD has the same blind spot. An iteration that only multiplies computes the polar factor by matrix products alone, and each singular value follows its own scalar map; nothing in that map can see the determinant, so Newton–Schulz converges to the mirrored polar factor whenever the SVD would have returned it. A rotation from the inverse-free route needs the determinant checked afterwards and the smallest direction flipped by some other means — which needs the smallest singular vector, the one object the iteration was chosen to avoid computing.

And the question is exactly the one rank is a decision raises about a different small singular value. There, a σ near rounding is a decision about whether a direction exists; here, a σ₃ near the noise is a decision about which way a direction points. In both cases the arithmetic returns a definite answer, and the definiteness is not evidence.

Where the polar factor goes from here #

Weights and the Wahba form. Attitude determination weights each sighting by its reliability, which changes M but not the flip rule. Whether weighting a few precise points heavily behaves like increasing m in the second figure — raising the effective thickness — or like a smaller point set, with its departure from the √m pattern, is measurable with the same computation.

Dimensions above three. The rule generalises unchanged — flip the direction of the smallest singular value — and the probability that the polar factor is improper should depend on the smallest singular value alone, but the number of directions a noisy estimate can be thin in grows with n. Registration in higher dimensions, and the orthogonal factor of a latent-space alignment, are where a mirror would go unnoticed longest, since there is no picture to look at.

A sign-aware iteration. The inverse-free iteration’s blind spot is a scalar map with fixed points at ±1. An iteration that converged to the rotation directly — carrying a determinant, or deflating the smallest direction as it goes — would keep the multiply-only property and remove the need for an SVD at the end. Whether one exists with Newton–Schulz’s √3 basin, or whether the determinant constraint necessarily costs a factorisation, is the question this measurement leaves for the next.