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²  =  ‖Σ − PᵀWQ‖_F²

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.

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₂‖_F  =  3.10,   on matrices whose own norm is 2.83

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.

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 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.

‖UᵀU − I‖_F      the orthogonality of the factor
‖A − UH‖_F/‖A‖   the reconstruction

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 σ_k, 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.