The nearest orthogonal matrix

A rotation matrix in a robot’s state estimator drifts. A basis in a subspace iteration stops being orthonormal. A transformation fitted to noisy point pairs comes out not quite orthogonal. In every case somebody has a matrix that is nearly orthogonal and wants the orthogonal matrix nearest to it, and in every case the answer written down is QR.

It is the wrong answer, and the interesting part is not by how much.

What the nearest orthogonal matrix is #

Write A = PΣQᵀ, its singular value decomposition. Then

A  =  (PQᵀ)(QΣQᵀ)  =  U H

with U = PQᵀ orthogonal and H = QΣQᵀ symmetric positive semidefinite. That is the polar decomposition, and U is the nearest orthogonal matrix to A in the Frobenius norm, in the spectral norm, and in every unitarily invariant norm — one matrix minimising all of them at once, which is a stronger statement than it looks.

The proof is one line once the SVD is in hand: for orthogonal W,

$$\|A - W\|_F^2 = \|\Sigma - P^{\mathsf T} W Q\|_F^2$$

and PᵀWQ is orthogonal, so the closest it can come to a diagonal of positive numbers is the identity, which is W = PQᵀ. Nothing to choose.

Fifty-five times is the largest ratio anywhere in this essay and it is on the easiest matrix in it, which is the first sign that one of these four candidates is not being measured the way the other three are:

Two of the four candidates have moved a long way and two have barely moved, and one of the four has moved in the opposite direction to the rest — which is not what the minimisation guarantees and is the correction this section ends on.

One decade further and the random baseline overtakes the cheap method, which is worth watching for because it is the point at which the whole question stops paying:

Correction to that caption: not every candidate moves further away. Polar and sign-fixed QR do — 0.08 to 2.52 and 0.09 to 2.87 — but raw QR goes the other way, 4.4225 down to 3.1329, and the random baseline goes 3.1110 down to 2.7435. Only the two methods that are actually solving the minimisation get worse as the matrix gets harder. The two that are not solving it converge on the same place from above.

Sign-fixed QR is within 15% of optimal at every κ: 1.099, 1.146, 1.140, 1.133, 1.138. Five conditionings, and the cheap route’s penalty is a constant to two figures. That is the practical result of the essay and it is stronger stated across the slider than at one point.

The sign fix matters most where the matrix needs the least correcting. Raw QR is 54.7× out at κ = 1.05 and 1.24× out at κ = 10⁴, because a flipped column costs a fixed amount and the amount being measured against is 0.08 at one end and 2.5 at the other. So the failure to fix signs is invisible on the hard problems and catastrophic on the easy ones, which is the wrong way round for anybody testing an implementation.

And the problem becomes vacuous as κ grows. At κ = 1.05 the optimal answer is 38 times better than the best of two hundred random orthogonal matrices; at κ = 10⁴ it is better by 8.9% — 2.5188 against 2.7435. A badly stretched matrix has no nearby orthogonal matrix, so nearest stops being a useful word, and the two hundred random draws are there to say so rather than as decoration. Raw QR, meanwhile, is beaten by those random draws at every κ on the slider.

The measurement, including the two hundred controls #

The strip beneath the bars is two hundred randomly drawn orthogonal matrices, and it is there because nearest is a claim about all of them.

The polar factor is nearer at every conditioning, by between 9 and 17 per cent. That is a real gap and it is a small one, and it would not on its own be worth an essay.

The gap that is not small #

Permute the columns of A and ask the same question again.

polar(AP)  =  polar(A) P     to 2.3·10⁻¹⁵
QR:   ‖Q₁P − Q₂‖  =  3.10,     on matrices whose own norm is 2.83, both in the Frobenius norm

QR’s answer to orthogonalise this depends on the order the columns were written in, and the polar factor’s does not. The two QR answers are not perturbations of each other; they are unrelated matrices, differing by more than their own size.

The reason is structural rather than numerical, and it is visible in the algorithms. Gram–Schmidt leaves the first column alone and projects the last against everything before it. Householder introduces zeros one column at a time. Both are sequential in the columns, so they treat the first one differently from the last, and permuting the columns changes which column is which.

The polar factor is defined by a minimisation over the whole matrix at once. There is no order in it to depend on.

The previous phase on this site established a test for exactly this situation: a number that moves under something the problem does not move under is a measure of the description rather than of the problem. Permuting columns does not change the column space, the singular values, or which orthogonal matrix is nearest. QR’s answer moves anyway, so QR is answering a question about the ordering.

The gap does not close where the method is used #

The table starts at κ = 1.05, and every case this essay opens with is much nearer orthogonal than that: a rotation that has drifted, a basis that has stopped being orthonormal, a transformation fitted to slightly noisy pairs. All of them sit at κ − 1 of a thousandth or less, and the natural expectation is that the two answers converge there and QR becomes harmless.

They do not converge. Taking κ down to 1.001:

at κ = 1.001 the polar factor is 1.689·10⁻³ from A and QR’s Q is 1.861·10⁻³ — an excess of 10.2 per cent. At 1.01, 1.675·10⁻² against 1.846·10⁻², again 10.2 per cent. At 1.05, 10.2 per cent. At 1.5, 10.7. At 50, 11.8. At 1000, 9.8.

The excess is ten per cent at every conditioning, over five orders of magnitude in κ − 1. It is not a property of hard matrices and it does not fade in the regime the method is actually used in.

And the ambiguity is the size of the correction #

The permutation gap is the sharper measurement once it is put in the right units.

The essay compares it to the matrix’s own norm — 3.10 against 2.83 — which says it is large. The comparison that decides whether it matters is against the correction being computed, because that is what the answer is for: A is nearly orthogonal, and what is wanted is the small change that makes it orthogonal exactly.

Measured, ‖Q(A)P − Q(AP)‖ divided by ‖A − Q(A)‖: 0.84 at κ = 1.001, 0.84 at 1.01, 0.86 at 1.05, 0.97 at 1.5. The ratio is the same at every conditioning, and it is close to one.

So the two QR answers differ from each other by very nearly as much as either differs from A. The arbitrary part of the correction is the size of the correction. A drifted rotation handed to QR gets moved by some amount, and permuting its columns first moves it by a comparable amount in an unrelated direction — which is not a small perturbation of an answer, it is an answer with no more signal than noise in it.

That is what turns the ordering dependence from a curiosity into the finding. Ten per cent of extra distance is a defensible price for a cheaper factorisation, and if that were the whole story the sensible conclusion would be to use QR and not worry. An answer whose arbitrary half is as large as its useful half is a different object: it is not a worse correction, it is not reliably a correction at all.

And the ratio being constant is what says the mechanism is structural rather than a badly conditioned corner. If the ambiguity came from the matrix being hard, it would grow with κ against a correction that also grows, and the ratio would move. It does not move: the ordering dependence is proportional to how far A is from orthogonal, which is exactly what a sequential algorithm applied to a nearly orthogonal matrix must produce — each column is orthogonalised against its predecessors, the predecessors are different after a permutation, and the difference is first order in the drift.

A sign convention worth seven times the gap #

While measuring the above, something went wrong that is worth publishing rather than quietly fixing.

qrHouseholder returns an R whose diagonal entries may be negative, which negates the corresponding columns of Q. Both factorisations are perfectly valid — ‖A − QR‖/‖A‖ is 10⁻¹⁶ either way — and the convention is invisible in every quantity anybody checks.

It is not invisible in ‖A − Q‖. On the κ = 1.5 matrix, with five of eight columns flipped:

polar factor:        0.586
QR, signs fixed:     0.650
QR, as returned:     4.130
best of 200 random:  2.956

The unfixed Q is further from A than a randomly chosen orthogonal matrix. A figure drawn from it would have made QR look ridiculous, and the reason would have had nothing to do with QR.

Nothing in this fleet checks a sign convention. The factorisation is correct, the residual is at rounding, the columns are orthonormal, and every assertion the site has ever written about a QR passes. It is the same shape as the defect the last phase found in a marker’s radius: a property nobody thought to state, on an object every gate approves of.

What QR is for, which is something else #

None of this makes QR the wrong tool. It makes it the right tool for a different job, and the two jobs are worth separating because they are usually conflated under the word orthogonalise.

QR answers a question about a nested sequence of subspaces. Its Q has the property that the first k columns span the same space as the first k columns of A, for every k. That is a strong and useful property — it is why QR is the right factorisation for least squares, where the columns arrive in a meaningful order and the triangular R is the object being solved with — and the polar factor does not have it at all.

The polar factor answers a question about the whole matrix at once. Nearest, in every unitarily invariant norm, with no reference to any ordering. It has no triangular partner and it solves nothing.

So the rule is: if the columns mean something in order, QR; if the matrix means something as a whole, polar. A least-squares design matrix is the first case. A rotation that has drifted is the second. Confusing them is what produces an answer that changes when somebody reorders the inputs.

There is a third case worth naming, because it is the one where QR is reached for and is genuinely inadequate rather than merely different: a basis for a subspace that is going to be iterated on. Subspace iteration re-orthogonalises its basis every step, the columns have no meaningful order, and QR’s answer therefore depends on an arbitrary convention that changes as the basis rotates. Everything converges anyway, and the intermediate quantities are noisier than they need to be.

Where the answer is actually wanted #

The polar factor is the answer to a question three fields ask independently and each named separately.

Orthogonal Procrustes. Given two point sets, find the rotation that best takes one onto the other: minimise ‖A − BQ‖_F over orthogonal Q. The answer is the polar factor of BᵀA. In structural biology this is the Kabsch algorithm and it aligns protein structures; in spacecraft attitude determination it is Wahba’s problem; in computer vision it is the rotation step of iterative closest point. Same formula, three names, three literatures.

Re-orthogonalising a drifted rotation. A quaternion or a rotation matrix integrated forward in time stops being orthogonal — the same drift a sequence accumulates one field over. The nearest orthogonal matrix is the right correction, and it is the one that moves the estimate least — which matters, because the correction is being applied continuously and any bias in it accumulates.

Density matrix purification and the sign function. In electronic structure calculations the idempotent projector nearest a given symmetric matrix is built from the same iteration, and the subject is the same one.

Two routes to the same U #

The polar factor above is computed through an SVD, which is honest and expensive: an SVD costs more than the eigendecomposition it is often being used to avoid. It is worth checking the answer a second way before spending an essay on cheaper routes.

The second route is the SVD’s own definition run backwards. U = PQᵀ is orthogonal, so UᵀU = I, and H = UᵀA must come out symmetric positive definite — two properties neither of which was imposed. Both hold to 10⁻¹⁵ on every matrix measured here, and the symmetry of UᵀA is a genuinely independent check: nothing in the construction of PQᵀ makes QΣQᵀ symmetric except the algebra being right.

A third route exists and is the subject of the next essay: iterate. The interesting thing about the iterations is that they never mention Σ, P or Q at all, so where they agree with PQᵀ the agreement is between two computations with almost nothing in common.

The residual, printed #

The site’s rule applies here as everywhere: no decomposition drawn without its residual. For the polar decomposition there are two.

Both are at 10⁻¹⁵ throughout the measurements above. That matters because the polar factor is computed here through an SVD, and an SVD is an expensive thing to insist on — the next essay is entirely about two ways of getting U without one, and the residuals are how those two ways are held to the same standard as this one.

What H is, which nobody asks about #

The polar decomposition has a second factor and it is usually thrown away. It should not be, because it is the more interpretable of the two.

H = QΣQᵀ is symmetric positive semidefinite with exactly the singular values of A as its eigenvalues. So A = UH says: a linear map is a stretch along orthogonal axes, followed by a rotation. That is the geometric statement the singular value decomposition is usually explained with, and the polar form says it with two objects instead of three.

In continuum mechanics the same decomposition of the deformation gradient is the right stretch tensor and the rotation, and separating them is what makes a material law frame-indifferent — the stress must depend on H and not on U, because U is where the observer’s choice of orientation went. An entire field’s constitutive theory rests on the factorisation this essay is about, reached independently and for a reason that has nothing to do with numerical analysis: U carries the part of the map that a change of observer changes, and H carries the part that it does not.

That is the same invariance argument as the column-permutation one above, and it is worth noticing that it comes from mechanics rather than from computation. A quantity that moves under something the physics does not move under cannot be in a material law; a quantity that moves under something the question does not move under cannot be the answer to it.

Rotations, reflections, and the constraint nobody states #

There is a constraint the minimisation above does not carry and most applications need, and it is the source of a well-known bug.

The orthogonal matrices come in two components: det = +1, the rotations, and det = −1, the reflections. The polar factor minimises ‖A − Q‖ over all of them, so on a data set that has been reflected — by a sign convention, a left-handed coordinate frame, or noise on a nearly degenerate configuration — it returns the reflection, which is the correct answer to the question asked and the wrong answer to the question meant. A molecular structure comes back mirrored; a robot’s estimated attitude flips.

The repair is in every implementation of the Kabsch algorithm and is worth writing out because it is not obvious that it is optimal. Compute the SVD of BᵀA = PΣQᵀ, and if det(PQᵀ) < 0, replace it by P·diag(1, …, 1, −1)·Qᵀ — that is, flip the sign of the singular vector belonging to the smallest singular value. The result is the nearest matrix with determinant +1, and the reason it is the smallest one to flip is the same argument as the original minimisation: the cost of flipping the kth direction is proportional to σₖ, so flipping the cheapest is optimal.

That is a genuinely satisfying result and it is also a warning. The unconstrained answer and the constrained one differ by a reflection, which is not a small perturbation of anything. A code that uses the polar factor for a rotation and does not check the determinant is correct on every input where the data is not nearly degenerate, and catastrophically wrong on the ones where it is — which is the shape of every intermittent bug in this fleet.

And the determinant makes an appearance the phase’s first essay would not have predicted: it is the right test here because the input is orthogonal, so its determinant is exactly ±1 and the comparison is against zero with no threshold in it. The scalar this phase opened by arguing against is the correct diagnostic on the one class of matrices where it carries no units at all.

What is worth carrying #

QR does not answer the question it is usually reached for. The nearest orthogonal matrix is the polar factor, nearer by 9 to 17 per cent — and the percentage is not the reason to use it.

The reason is that the answer does not depend on the column order. polar(AP) = polar(A)P to fourteen digits; QR’s answer under the same permutation is a different matrix by more than its own size. If a question does not care about the order of the columns, an answer that does is answering a different question.

And a sign convention is worth seven times the gap it hides. The raw Householder Q is further from A than a random orthogonal matrix, both factorisations reconstruct A perfectly, and no assertion on this site would have caught it.

The next essay computes the same U without an SVD, twice, and finds that one of the two methods converges to an orthogonal matrix that is not the answer and says nothing about it: an iteration that only multiplies.