Maths - Reorthogonalising a matrix

We may start with orthogonal matrices, but when we begin to do operations on them small rounding errors can build up, causing the matrices to be slightly de-orthogonalised. Therefore, after many such operations, we may need a method to reorthonalise them.

where we doing a lot of calculations and errors can build up, for

An example of where reorthogonalising a matrix may be important is when we are recalculating at every frame in an animation, as described here.

Unfortunately I have not found an easy formula for doing this, here are the methods I am aware of:

None of these methods are ideal in that they either seem to have inaccuracies or they are very complicated, I would welcome any suggestions for better methods or improvements.

Recalculating Basis Vectors

The idea here is to remove the redundancy involved in a rotation matrix and then reconstruct the missing bits so that the matrix is perfectly orthogonal.

This is the simplest approach, but it is not ideal in that the redundant information might have helped us cancel out the rounding errors, instead this approach may degrade the accuracy of the result but the matrix will be orthogonal.

If we consider a rotation matrix as a set of basis vectors (as discussed on this page) then we need to make sure these basis vectors are unit length, divide by √(x² + y² + z²), and also mutually perpendicular. We could do that by deleting the third basis vector and recalculating it from the cross product of the other two. So far we have not made sure that the first two basis vectors are perpendicular with each other, so we the need to delete on of them and then recalculate from the other two.

Another approach could be to convert the matrix to a quaternion and then convert back to a matrix.

Either way its messy and we could be basing our results on the least accurate elements of the matrix.

Deriving Correction Matrix

The following derivation is evolved from this discussion Colin Plumb on newsgroup: sci.math

We want a correction matrix [C] which when multiplied with the original matrix [M] to give an orthogonal matrix [O] such that:

[O] = [C]*[M]

and since [O] is orthogonal then,

[O]-1 = [O]T

combining these equations gives:

[C]*[M] * ([C]*[M])T = [I]

where [I] = identity matrix

as explained here, ([A] * [B])T = [B]T * [A]T so substituting for [C] and [M] gives,

[C]*[M] * [M]T*[C]T = [I]

post multiplying both sides by ([C]T)-1 gives:

[C]*[M] * [M]T = ([C]T)-1

pre multiplying both sides by [C]-1 gives:

[M] * [M]T = [C]-1 ([C]T)-1

If we constrain [C] to be symmetric then [C]T = [C]

[M] * [M]T = [C]-2

[C] = ([M] * [M]T)-(1/2)

So using the taylor series here and using the first two terms we get:

x^(-1/2) = (3-x)/2

So substituting in above gives:

[C] = ([M] * [M]T)-(1/2)

= (3-([M] * [M]T))/2

since we can't subtract a matrix from a scalar I think it should really be written:

[C] = (3[I]-([M] * [M]T))/2

At least this seems to give appropriate results in the examples below:

There are some things about this derivation that make me uneasy, for instance constraining [C] to be symmetric seems to come out of nowhere. Vassilis has raised some issues with this method here, can anyone help me resolve them?

SVD Method

Hauke kindly suggested this method.

It is based on factoring the matrix into 3 parts:

[M] = [U][D][V]t

where:

More information about SVD on this page.

This has separated the matrix our into a scaling part, the diagonal matrix [D] and the rotational part [U] and [V].

We want no (unity) scaling so it is tempting to set all the terms in the diagonal matrix to one:

[D] = [I] = identity matrix

which would give:

[M] = [U][V]t

However, because of the way that the SVD algorithm calculates [D], [U] and [V] it is possible that the terms of [D] might be -1 as well as +1. We therefore need to set the terms of [D] to +1 if the term was positive and -1 if the term was negative.

Alternatively we could recalculate the diagonal terms to give a +1 determinant.

D=diag([1 1 det(U*V^T)]).

Which means that our orthogonal matrix is:

[M] = U*diag([1 1 det(U*V^T])*V^T.

Examples

Orthogonal matrix

First test that this will preserve the orthogonal matrix, applying:

= (3[I]-([M] * [M]T))/2

Since the matrix is orthogonal then [M] * [M]T = [I] so,

[C] = (3[I]-(
1 0
0 1
))/2

 

[C] =
(3- 1)/2 0
0 (3- 1)/2

 

[C] =
1 0
0 1

so as required this wont alter the identity matrix as required:

Non-orthogonal matrix

now try this on a matrix that is not already orthogonal because it has a term 'error' added to m01, applying:

[C]= (3-([M] * [M]T))/2

[C]= (3-(
1 error
0 1
1 0
error 1
))/2

 

[C]=
(3-1+error*error))/2 -error/2
-error/2 (3- 1)/2

 

[C]=
1 + 3*error*error/2 -error/2
-error/2 1

so if we multiply this be the original matrix we get:

[O]=
1 + 3*error*error/2 -error/2
-error/2 1
1 error
0 1

 

[O]=
1 + 3*error*error/2 error + 3*error*error*error/2 - error/2
-error/2 1-error*error/2

 

[O]=
1 + 3*error*error/2 (error + 3*error*error*error)/2
-error/2 1-error*error/2

Is this more orthogonal than the original matrix, I can't really tell, it looks as though it might be?


metadata block
see also:
Correspondence about this page

Book Shop - Further reading.

Where I can, I have put links to Amazon for books that are relevant to the subject, click on the appropriate country flag to get more details of the book or to buy it from them.

cover Mathematics for 3D game Programming - Includes introduction to Vectors, Matrices, Transforms and Trigonometry. (But no euler angles or quaternions). Also includes ray tracing and some linear & rotational physics also collision detection (but not collision response).

Terminology and Notation

Specific to this page here:

 

This site may have errors. Don't use for critical systems.

Copyright (c) 1998-2016 Martin John Baker - All rights reserved - privacy policy.