From: kovarik@mcmail.cis.McMaster.CA (Zdislav V. Kovarik)
Subject: Re: renormalizing rotation matrices
Date: 25 Aug 2000 13:43:12 -0400
Newsgroups: sci.math.num-analysis
Summary: [missing]

In article <8o64vm$sde$1@nnrp1.deja.com>,  <mike_hayford@my-deja.com> wrote:
:I am working with the implementation of a CAD system where imported data
:sometimes has rotation matrices that are not orthonormal, either because
:of truncation of data or roundoff in previous calculations.  I would
:like to calculate the nearest orthonormal matrix when I detect problems.
:
:I can think of a number of adhoc approaches such as running the matrix
:thru a Gram-Schmidt orthogonalizer or doing an eignevalue decomposition
:and then setting the eignevalues to unity and reforming.  These seem
:to work on the problems I've tried so far but I don't know if or when
:they'll fail.  I'd be interested in pointers to any rigorously thought
:through approaches to the problem.

 I don't have the reference with me right now but it is an old result:
Given a matrix M which is close to orthogonal, the closest to M orthogonal
matrix W (in Frobenius norm) is (M' means M transpose)

    W = M * (M'*M)^(-1/2)

 If you have Singular Value Decomposition in your software, this can be
accomplished by a variant of what you have done:

    Suppose M = U * S * V'  with U,V orthogonal and S diagonal with
positive diagonal entries (they usually come sorted from the largest),
then the same W is

   W = U * V'   (setting the singular values to 1)

Iterative improvement, provided that (I - M'*M) has norm less than 1 (if M
comes as a small perturbation of an originally orthogonal matrix, this
should be the case):

   Set M(0)=M,  R(0) = I - M'*M

   given M(k), R(k), calculate M(k+1) = M(k)*(I + (1/2)*R(k))
                     and of course R(k+1) = I - M(k)' *M(k)

This iteration converges quadratically to W if done with full accuracy.

Even with small intermediate errors made at every step, the distance from
M(k) to an orthogonal matrix closest to M(k) is not more than

  r(k) / (1 + sqrt(1 - r(k))  where r(k) is the norm of R(k).

Hope it helps, ZVK(Slavek).