Thread Subject:

Procrustes Analysis without Reflection

Hi all,

I'm working on a problem that requires the procrustes alignment of two
sets of data, but it is important that the proposed transformation
does not allow reflection. The MatLab procrustes function allows this
and I was wondering if anyone knows of a freely available m-file that
will perform Procrustes analysis without allowing reflection.

Thanks in advance,

Ben.

rizwindu <bgeorge@gmail.com> wrote in message <435f98ff-
fb02-495f-9600-16fbad0846ac@w7g2000hsa.googlegroups.com>...
> Hi all,
>
> I'm working on a problem that requires the procrustes alignment of two
> sets of data, but it is important that the proposed transformation
> does not allow reflection. The MatLab procrustes function allows this
> and I was wondering if anyone knows of a freely available m-file that
> will perform Procrustes analysis without allowing reflection.
>
> Thanks in advance,
>
> Ben.
------------
  This must be a hasty reply because I have a doctor's appointment in just a
few minutes. Later in the day I can send a more complete response.

  If X and Y are the two n x m data sets in question where n is the number of
data items and m the number of components to the data, then one step of
the procrustes procedure is to do a singular value decomposition, 'svd', on
the m x m matrix X'*Y, getting

 X'*Y = U*S*V'

where U and V are unitary matrices and S is diagonal with the non-negative
singular values. The best fitting unitary transformation in the least squares
sense is given by Q = V*U'. However, if it turns out that Q involves a
reflection, that is, its determinant is -1, then what you need to do is reverse
the sign in the rightmost column of either U or V, corresponding to the
smallest singular value. Then do Q = V*U' with the modified U or V, and you
get a unitary transformation without reflection which is the best fit among
such non-reflective transformations. The translation and scaling steps
remain the same.

  It should be said that if you would have had a reasonably close fit using a
reflected transformation and the smallest singular value is fairly large, then
you will get a comparatively poor fit with the best non-reflected
transformation. It would be analogous to trying to fit the corresponding
vertices of a fat triangle onto its mirror image without turning it over.

  More later. I must go now. My apologies.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g09nco$6f1$1@fred.mathworks.com>...
> This must be a hasty reply because I have a doctor's appointment in just a
> few minutes. Later in the day I can send a more complete response.
>
> If X and Y are the two n x m data sets in question where n is the number
of
> data items and m the number of components to the data, then one step of
> the procrustes procedure is to do a singular value decomposition, 'svd', on
> the m x m matrix X'*Y, getting
>
> X'*Y = U*S*V'
>
> where U and V are unitary matrices and S is diagonal with the non-negative
> singular values. The best fitting unitary transformation in the least squares
> sense is given by Q = V*U'. However, if it turns out that Q involves a
> reflection, that is, its determinant is -1, then what you need to do is
reverse
> the sign in the rightmost column of either U or V, corresponding to the
> smallest singular value. Then do Q = V*U' with the modified U or V, and
you
> get a unitary transformation without reflection which is the best fit among
> such non-reflective transformations. The translation and scaling steps
> remain the same.
>
> It should be said that if you would have had a reasonably close fit using a
> reflected transformation and the smallest singular value is fairly large, then
> you will get a comparatively poor fit with the best non-reflected
> transformation. It would be analogous to trying to fit the corresponding
> vertices of a fat triangle onto its mirror image without turning it over.
>
> More later. I must go now. My apologies.
>
> Roger Stafford
------------
  I'm back from the doctor and will try to explain what I did earlier in greater
detail. You have two n x m arrays of data, X and Y. You desire to translate,
rotate, and rescale the data in X so as to best match that in Y in the least
squares sense. It is easy to show that the best translation is to have the
mean value of the X data equal to that of Y, so we now assume that both
means are located at the origin. This allows us to then rotate about the
origin.

  The rotation is more difficult. First we seek a unitary matrix Q which will
transform X to best match Y. The squares-sum quantity we wish to minimize
can be expressed as:

 ||X*Q - Y||^2 = trace((X*Q-Y)'*(X*Q-Y)) =
 trace(Q'*X'*X*Q) + trace(Y'*Y) - trace((X*Q)'*Y) - trace(Y'*(X*Q)) =
 trace(Q*Q'*X'*X) + trace(Y'*Y) - trace(Q'*X'*Y) - trace(Y'*X*Q) =
 trace(X'*X) + trace(Y'Y) - 2*trace(Y'*X*Q)

Hence we wish to select Q so as to maximize trace(Y'*X*Q). The matrix Q is
understood to be unitary - that is, its transpose is its inverse. The
determinant of Q can be either +1 or -1. If it is +1, it is a valid orthogonal
rotation with no reflection involved. If, for the moment, we are unconcerned
whether Q is a rotation, the best Q can be determined using matlab's singular
value decomposition function, 'svd':

 Y'*X = U*S*V'

where U and V are unitary and S has the singular values in descending order
along the diagonal. This gives

 trace(Y'*X*Q) = trace(U*S*V'*Q) = trace(S*V'*Q*U)

The quantity V'*Q*U is also unitary and we get the maximum value for this
trace if it is set equal to I, the m x m identity matrix. Thus

 trace(S*I) = trace(S) = sum of the singular values,

is the maximum value for the above, and it is obtained by setting

 V'*Q*U = I, or
 Q = V*U'.

  However, V and U may be such that det(Q) = -1 and in accordance with your
stated condition, this cannot be used. The remedy is, as I have said earlier,
to reverse the sign on the rightmost column of either U or V where the
smallest singular value is located. When that is done, we have

 Y'*X = U2*S2*V2'

where U2 or V2 has been so altered and where S2 is S except that the last
diagonal value is the negative of that in S. In other words we have a slightly
altered form of singular value decomposition in which the smallest "singular
value" is made to be negative. In that case, the above trace which we were
maximizing becomes

 trace(S2*I) = trace(S2) =
  sum of first m-1 singular values minus the last, smallest one

with Q = V2*U2', which will now be a valid rotation with determinant +1.

  I confess I don't have a rigorous proof that this last is actually the maximum
possible value for m > 2, but I ran some rather exhaustive tests which
verified that it is true for m = 3. My intuition tells me it is also true for all m,
but for now it will have to remain as strongly-held intuition.

  The last step is rescaling. What scale factor, f, can we multiply X*Q by to
achieve a final minimization? If we write

 trace((f*X*Q-Y)'*(f*X*Q-Y)) =
 f^2*trace(X'*X) + trace(Y'Y) - 2*f*trace(Y'*X*Q)

it is evident that we must have

 f = trace(S)/trace(X'*X) or trace(S2)/trace(X'*X)

for the minimum value and therefore the closest least squares match. I
believe I stated earlier that rescaling remained the same, but as you see here,
it is actually somewhat changed.

Roger Stafford

I am trying to understand exactly what is done in Matlab's
PROCRUSTES. From your reply, it seems that given two sets
of points, it first assumes that the center of rotation is
at the mean of each set, which reduces the problem to just
finding the rotation matrix using SVD. I am dealing with a
more general problem in which I cannot assume that the sets
of points are symmetrical around the mean. In that case,
the mean assumption would lead to grossly inaccurate
estimates of the rotation matrix. In the general problem,
the center-of-rotation vector just adds N unknowns to an
N-dim shape. In principle given a succession of rotating
shapes plus random noise (i.e. more than N samples), it
should be possible to simultaneously estimate the
center-of-rotation and the rotation matrix. Is there any
way to do this simply in Matlab using a variant of
PROCRUSTES, or do I need to work this all out by myself?
Thanks
QZ


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <g0anam$kor$1@fred.mathworks.com>...
> "Roger Stafford"
<ellieandrogerxyzzy@mindspring.com.invalid> wrote in
> message <g09nco$6f1$1@fred.mathworks.com>...
> > This must be a hasty reply because I have a doctor's
appointment in just a
> > few minutes. Later in the day I can send a more
complete response.
> >
> > If X and Y are the two n x m data sets in question
where n is the number
> of
> > data items and m the number of components to the data,
then one step of
> > the procrustes procedure is to do a singular value
decomposition, 'svd', on
> > the m x m matrix X'*Y, getting
> >
> > X'*Y = U*S*V'
> >
> > where U and V are unitary matrices and S is diagonal
with the non-negative
> > singular values. The best fitting unitary
transformation in the least squares
> > sense is given by Q = V*U'. However, if it turns out
that Q involves a
> > reflection, that is, its determinant is -1, then what
you need to do is
> reverse
> > the sign in the rightmost column of either U or V,
corresponding to the
> > smallest singular value. Then do Q = V*U' with the
modified U or V, and
> you
> > get a unitary transformation without reflection which is
the best fit among
> > such non-reflective transformations. The translation
and scaling steps
> > remain the same.
> >
> > It should be said that if you would have had a
reasonably close fit using a
> > reflected transformation and the smallest singular value
is fairly large, then
> > you will get a comparatively poor fit with the best
non-reflected
> > transformation. It would be analogous to trying to fit
the corresponding
> > vertices of a fat triangle onto its mirror image without
turning it over.
> >
> > More later. I must go now. My apologies.
> >
> > Roger Stafford
> ------------
> I'm back from the doctor and will try to explain what I
did earlier in greater
> detail. You have two n x m arrays of data, X and Y. You
desire to translate,
> rotate, and rescale the data in X so as to best match that
in Y in the least
> squares sense. It is easy to show that the best
translation is to have the
> mean value of the X data equal to that of Y, so we now
assume that both
> means are located at the origin. This allows us to then
rotate about the
> origin.
>
> The rotation is more difficult. First we seek a unitary
matrix Q which will
> transform X to best match Y. The squares-sum quantity we
wish to minimize
> can be expressed as:
>
> ||X*Q - Y||^2 = trace((X*Q-Y)'*(X*Q-Y)) =
> trace(Q'*X'*X*Q) + trace(Y'*Y) - trace((X*Q)'*Y) -
trace(Y'*(X*Q)) =
> trace(Q*Q'*X'*X) + trace(Y'*Y) - trace(Q'*X'*Y) -
trace(Y'*X*Q) =
> trace(X'*X) + trace(Y'Y) - 2*trace(Y'*X*Q)
>
> Hence we wish to select Q so as to maximize trace(Y'*X*Q).
 The matrix Q is
> understood to be unitary - that is, its transpose is its
inverse. The
> determinant of Q can be either +1 or -1. If it is +1, it
is a valid orthogonal
> rotation with no reflection involved. If, for the moment,
we are unconcerned
> whether Q is a rotation, the best Q can be determined
using matlab's singular
> value decomposition function, 'svd':
>
> Y'*X = U*S*V'
>
> where U and V are unitary and S has the singular values in
descending order
> along the diagonal. This gives
>
> trace(Y'*X*Q) = trace(U*S*V'*Q) = trace(S*V'*Q*U)
>
> The quantity V'*Q*U is also unitary and we get the maximum
value for this
> trace if it is set equal to I, the m x m identity matrix.
 Thus
>
> trace(S*I) = trace(S) = sum of the singular values,
>
> is the maximum value for the above, and it is obtained by
setting
>
> V'*Q*U = I, or
> Q = V*U'.
>
> However, V and U may be such that det(Q) = -1 and in
accordance with your
> stated condition, this cannot be used. The remedy is, as
I have said earlier,
> to reverse the sign on the rightmost column of either U or
V where the
> smallest singular value is located. When that is done, we
have
>
> Y'*X = U2*S2*V2'
>
> where U2 or V2 has been so altered and where S2 is S
except that the last
> diagonal value is the negative of that in S. In other
words we have a slightly
> altered form of singular value decomposition in which the
smallest "singular
> value" is made to be negative. In that case, the above
trace which we were
> maximizing becomes
>
> trace(S2*I) = trace(S2) =
> sum of first m-1 singular values minus the last,
smallest one
>
> with Q = V2*U2', which will now be a valid rotation with
determinant +1.
>
> I confess I don't have a rigorous proof that this last
is actually the maximum
> possible value for m > 2, but I ran some rather exhaustive
tests which
> verified that it is true for m = 3. My intuition tells me
it is also true for all m,
> but for now it will have to remain as strongly-held intuition.
>
> The last step is rescaling. What scale factor, f, can
we multiply X*Q by to
> achieve a final minimization? If we write
>
> trace((f*X*Q-Y)'*(f*X*Q-Y)) =
> f^2*trace(X'*X) + trace(Y'Y) - 2*f*trace(Y'*X*Q)
>
> it is evident that we must have
>
> f = trace(S)/trace(X'*X) or trace(S2)/trace(X'*X)
>
> for the minimum value and therefore the closest least
squares match. I
> believe I stated earlier that rescaling remained the same,
but as you see here,
> it is actually somewhat changed.
>
> Roger Stafford
>

"Qasim " <qzqzqzqz@gmail.com> wrote in message <g55807$io8
$1@fred.mathworks.com>...
> I am trying to understand exactly what is done in Matlab's
> PROCRUSTES. From your reply, it seems that given two sets
> of points, it first assumes that the center of rotation is
> at the mean of each set, which reduces the problem to just
> finding the rotation matrix using SVD. I am dealing with a
> more general problem in which I cannot assume that the sets
> of points are symmetrical around the mean. In that case,
> the mean assumption would lead to grossly inaccurate
> estimates of the rotation matrix. In the general problem,
> the center-of-rotation vector just adds N unknowns to an
> N-dim shape. In principle given a succession of rotating
> shapes plus random noise (i.e. more than N samples), it
> should be possible to simultaneously estimate the
> center-of-rotation and the rotation matrix. Is there any
> way to do this simply in Matlab using a variant of
> PROCRUSTES, or do I need to work this all out by myself?
> Thanks
> QZ

  I'm afraid my wordage in this thread back in May 13 may have been
misleading, Qasim. Matlab's 'procrustes' does NOT assume that the mean
points of the two sets are coincident. Its first step is to translate one set in
such a way that these do become coincident.

  The original poster, Rizwindu, wished to align two sets of data in the same
way as 'procrustes' does, except that a possible reflection should be avoided.
It was not to provide a reflection even if that would furnish a better match. I
know of no way of forcing 'procrustes' to do that, so my endeavor in
answering this question was to show how one could achieve this goal without
using 'procrustes'.

  The three steps one would need to use would be 1) translate one of the sets
so that their mean points are coincident, 2) find the best fitting rotation of
one set about that common mean point so as to match the other set, and 3)
provide the best rescaling relative to this mean point. I went through 1) and
3) rather hastily because it is step 2) that is the most difficult. I explained
how step 2) could be modified in such a way as to prevent any reflection
occurring.

  As I outlined the procedure, for convenience I actually assumed that both
sets would have been translated so that their mean points lay at a common
origin, but this is very easy to do. You simply subtract each set's mean point
vector from each of its points so that its new mean then lies at the origin. It
is fairly easy to show that any least squares solution must involve a
translation that matches the two mean points followed by a rotation about
that common mean point.

  However the hard part is finding the correct rotation. If you have any
specific questions about that material I sent in May, please don't hesitate to
ask.

Roger Stafford

Hi Roger
Thanks for replying to my query. Consider this simple trial
program

u = normrnd(0,pi,[10,1]);
 v=sin(u)+2;
  X = [u,v];
  A = [-3.5 7.5];
  AA=repmat(A, 10,1);
  XX = X-AA;
 theta = pi*.2;
  S = [cos(theta) -sin(theta); sin(theta) cos(theta)];
   YY = XX*S;
  XXX = XX + normrnd(0,.3,[10 2]);
   Y=YY+AA+ normrnd(0,.3,[10 2])
   [d,Z,tr] = procrustes(XXX,Y);
  plot(XXX(:,1),XXX(:,2),'rx',...
    Y(:,1),Y(:,2),'b.',...
    Z(:,1),Z(:,2),'bx');
d
tr.b
tr.T
tr.c:


Essentially I define a shape X by 10 points that are the
sines of a randomly sampled interval. I then rotate it
around A which is explicitly very different from the mean
vector of the shape. The original shape X is distorted into
XXX with frame-specific random noise, as independently is
the rotated shape Y. When I ask procrustes to match Y and
XXX, it does quite a good job. It seems unlikely to me that
it is fitting a rotation matrix around the mean vector
(although I may just be dense), so I just want to understand
what method it is using. If you can help, that would be
great.

BTW people who need a reference for the method you were
describing in your initial reply, can check out "Challis,
J.H. (1995). A procedure for determining rigid body
transformation parameters. J. Biomechanics 28, 733-737."

Best wishes
QZ


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <g55kdg$omj$1@fred.mathworks.com>...
> "Qasim " <qzqzqzqz@gmail.com> wrote in message <g55807$io8
> $1@fred.mathworks.com>...
> > I am trying to understand exactly what is done in Matlab's
> > PROCRUSTES. From your reply, it seems that given two sets
> > of points, it first assumes that the center of rotation is
> > at the mean of each set, which reduces the problem to just
> > finding the rotation matrix using SVD. I am dealing with a
> > more general problem in which I cannot assume that the sets
> > of points are symmetrical around the mean. In that case,
> > the mean assumption would lead to grossly inaccurate
> > estimates of the rotation matrix. In the general problem,
> > the center-of-rotation vector just adds N unknowns to an
> > N-dim shape. In principle given a succession of rotating
> > shapes plus random noise (i.e. more than N samples), it
> > should be possible to simultaneously estimate the
> > center-of-rotation and the rotation matrix. Is there any
> > way to do this simply in Matlab using a variant of
> > PROCRUSTES, or do I need to work this all out by myself?
> > Thanks
> > QZ
>
> I'm afraid my wordage in this thread back in May 13 may
have been
> misleading, Qasim. Matlab's 'procrustes' does NOT assume
that the mean
> points of the two sets are coincident. Its first step is
to translate one set in
> such a way that these do become coincident.
>
> The original poster, Rizwindu, wished to align two sets
of data in the same
> way as 'procrustes' does, except that a possible
reflection should be avoided.
> It was not to provide a reflection even if that would
furnish a better match. I
> know of no way of forcing 'procrustes' to do that, so my
endeavor in
> answering this question was to show how one could achieve
this goal without
> using 'procrustes'.
>
> The three steps one would need to use would be 1)
translate one of the sets
> so that their mean points are coincident, 2) find the best
fitting rotation of
> one set about that common mean point so as to match the
other set, and 3)
> provide the best rescaling relative to this mean point. I
went through 1) and
> 3) rather hastily because it is step 2) that is the most
difficult. I explained
> how step 2) could be modified in such a way as to prevent
any reflection
> occurring.
>
> As I outlined the procedure, for convenience I actually
assumed that both
> sets would have been translated so that their mean points
lay at a common
> origin, but this is very easy to do. You simply subtract
each set's mean point
> vector from each of its points so that its new mean then
lies at the origin. It
> is fairly easy to show that any least squares solution
must involve a
> translation that matches the two mean points followed by a
rotation about
> that common mean point.
>
> However the hard part is finding the correct rotation.
If you have any
> specific questions about that material I sent in May,
please don't hesitate to
> ask.
>
> Roger Stafford
>
>

"Qasim " <qzqzqzqz@gmail.com> wrote in message <g55sct$shv
$1@fred.mathworks.com>...
> Hi Roger
> Thanks for replying to my query. Consider this simple trial
> program
>
> u = normrnd(0,pi,[10,1]);
> v=sin(u)+2;
> X = [u,v];
> A = [-3.5 7.5];
> AA=repmat(A, 10,1);
> XX = X-AA;
> theta = pi*.2;
> S = [cos(theta) -sin(theta); sin(theta) cos(theta)];
> YY = XX*S;
> XXX = XX + normrnd(0,.3,[10 2]);
> Y=YY+AA+ normrnd(0,.3,[10 2])
> [d,Z,tr] = procrustes(XXX,Y);
> plot(XXX(:,1),XXX(:,2),'rx',...
> Y(:,1),Y(:,2),'b.',...
> Z(:,1),Z(:,2),'bx');
> d
> tr.b
> tr.T
> tr.c:
>
>
> Essentially I define a shape X by 10 points that are the
> sines of a randomly sampled interval. I then rotate it
> around A which is explicitly very different from the mean
> vector of the shape. The original shape X is distorted into
> XXX with frame-specific random noise, as independently is
> the rotated shape Y. When I ask procrustes to match Y and
> XXX, it does quite a good job. It seems unlikely to me that
> it is fitting a rotation matrix around the mean vector
> (although I may just be dense), so I just want to understand
> what method it is using. If you can help, that would be
> great.
>
> BTW people who need a reference for the method you were
> describing in your initial reply, can check out "Challis,
> J.H. (1995). A procedure for determining rigid body
> transformation parameters. J. Biomechanics 28, 733-737."
>
> Best wishes
> QZ

  In two or three dimensions, any rigid-body rotation about a point A is the
mathematical equivalent of the same rotation about any other point B
combined with an appropriate translation. This means that your pure
rotation can be achieved by the same rotation about the mean vector of a set
combined with some translation. As it turns out, the least mean square
transformation must necessarily be a translation which moves one mean
point to the other one, combined with some rotation about this common
mean point. As I said, that is easily demonstrated. This is clearly what
'procrustes' is doing.

Roger Stafford

"Qasim " <qzqzqzqz@gmail.com> wrote in message
<g55sct$shv$1@fred.mathworks.com>...
> Hi Roger
> Thanks for replying to my query. Consider this simple trial
> program
>
> u = normrnd(0,pi,[10,1]);
> v=sin(u)+2;
> X = [u,v];
> A = [-3.5 7.5];
> AA=repmat(A, 10,1);
> XX = X-AA;
> theta = pi*.2;
> S = [cos(theta) -sin(theta); sin(theta) cos(theta)];
> YY = XX*S;
> XXX = XX + normrnd(0,.3,[10 2]);
> Y=YY+AA+ normrnd(0,.3,[10 2])
> [d,Z,tr] = procrustes(XXX,Y);
> plot(XXX(:,1),XXX(:,2),'rx',...
> Y(:,1),Y(:,2),'b.',...
> Z(:,1),Z(:,2),'bx');
> d
> tr.b
> tr.T
> tr.c:
>
>
> Essentially I define a shape X by 10 points that are the
> sines of a randomly sampled interval. I then rotate it
> around A which is explicitly very different from the mean
> vector of the shape. The original shape X is distorted into
> XXX with frame-specific random noise, as independently is
> the rotated shape Y. When I ask procrustes to match Y and
> XXX, it does quite a good job. It seems unlikely to me that
> it is fitting a rotation matrix around the mean vector
> (although I may just be dense), so I just want to understand
> what method it is using. If you can help, that would be
> great.
>
> BTW people who need a reference for the method you were
> describing in your initial reply, can check out "Challis,
> J.H. (1995). A procedure for determining rigid body
> transformation parameters. J. Biomechanics 28, 733-737."
>
> Best wishes
> QZ

Qasim,

There's an infinite number of points about which one can
rotate - that's why it's easiest to work with centred
configurations, and why PROCRUSTES centres the two
configurations as its first step.

Use 'edit procrustes' to see what it's doing (basically it
finds the rotation matrix using the SVD method Roger
described above)

To answer the question in your other thread - the
dissimilarity measure that PROCRUSTES returns is just the
sum of square errors but standardised for scale (i.e. the
solution is just the 'least-squares' solution)

Best wishes, S

Hi to everybody.

I have the same problem of finding the optimal rotation matrix, but in a complex space. Hence the rotation matrix is actually a unitary matrix of size 2x2 of complex entries. Can I use the same procedure of Mr Stafford with the transpose-conjugate operator instead of only transpose?
Concerning the problem of determinant +/-1 of this post (see bottom), I found a very nice solution at:

http://www.kwon3d.com/theory/jkinem/rotmat.html

In my problem with complex data the determinant is exp(i*phi), hence I don't know how the previous problem of the determinant maps in my case. Is still a problem, or the proposed solution Q = V*U' is always fine in my case?

Thanks
Paolo

> If X and Y are the two n x m data sets in question where n is the number of
> data items and m the number of components to the data, then one step of
> the procrustes procedure is to do a singular value decomposition, 'svd', on
> the m x m matrix X'*Y, getting
>
> X'*Y = U*S*V'
>
> where U and V are unitary matrices and S is diagonal with the non-negative
> singular values. The best fitting unitary transformation in the least squares
> sense is given by Q = V*U'. However, if it turns out that Q involves a
> reflection, that is, its determinant is -1, then what you need to do is reverse
> the sign in the rightmost column of either U or V, corresponding to the
> smallest singular value.

Hi to everybody.

I have the same problem of finding the optimal rotation matrix, but in a complex space. Hence the rotation matrix is actually a unitary matrix of size 2x2 of complex entries. Can I use the same procedure of Mr Stafford with the transpose-conjugate operator instead of only transpose?
Concerning the problem of determinant +/-1 of this post (see bottom), I found a very nice solution at:

http://www.kwon3d.com/theory/jkinem/rotmat.html

In my problem with complex data the determinant is exp(i*phi), hence I don't know how the previous problem of the determinant maps in my case. Is still a problem, or the proposed solution Q = V*U' is always fine in my case?

Thanks
Paolo

> If X and Y are the two n x m data sets in question where n is the number of
> data items and m the number of components to the data, then one step of
> the procrustes procedure is to do a singular value decomposition, 'svd', on
> the m x m matrix X'*Y, getting
>
> X'*Y = U*S*V'
>
> where U and V are unitary matrices and S is diagonal with the non-negative
> singular values. The best fitting unitary transformation in the least squares
> sense is given by Q = V*U'. However, if it turns out that Q involves a
> reflection, that is, its determinant is -1, then what you need to do is reverse
> the sign in the rightmost column of either U or V, corresponding to the
> smallest singular value.