Thread Subject:

GPA using the procrustes function

If one has to align multiple sets of data, is it possibly to iteratively run the function procrustes to match against the mean of the dataset until there is some kind of convergence ? Or is there a Generalized Procrustes Analysis function somewhere ?

A postscript to my query. Iterative pairwise procrustes to a common reference does work.

"Shyamal L" <lshyamalNOSPAM@yahoo.com> wrote in message <ia6550$11n$1@fred.mathworks.com>...
> If one has to align multiple sets of data, is it possibly to iteratively run the function procrustes to match against the mean of the dataset until there is some kind of convergence ? Or is there a Generalized Procrustes Analysis function somewhere ?
- - - - - - - -
  I'm not sure of your meaning where you say "iteratively run the function procrustes to match against the mean of the dataset". Are you proposing to run each of the members of a data set against their single mean? That wouldn't give meaningful results in my opinion. Also I don't know how "convergence" would have a bearing here.

  When you run 'procrustes' on a pair of data sets, it finds the translation, unitary transformation, and scaling from one to the other that possesses the least mean square differences. The unitary transformation can involve a mirror imaging if that reduces the mean square differences. However taking the mean of a multiple set of data might be an almost meaningless entity to match with its members. For example suppose you have just two different 2D data sets, one the result of a 180 degree rotation of the other. Their mean would be a single point located at the center of rotation. What would be the point of matching each of the two members with a single point?

  I presume you have something different in mind but I haven't managed to figure out what it might be. Can you give a clearer description of your question? In particular tell just what you mean by "align multiple sets of data". Align them with what - with one of their members? If so, which one?

Roger Stafford

I know i could use one of the sets of shapes (k objects with m points having n dimensions) as my reference for aligning the rest. Fortunately it seems my data was not too off from each other so I used the mean of the k objects for each of the m points but I see your point that this could be quite an ugly and useless target to match against in the more general case.

My convergence criterion was merely the sum of squares between each pair of points for the k(k-1)/2 items and to see if the reduction at each iteration went below a threshold.

function Aligned = GPA(X)
[m,n,k] = size(X);
% assume k > 1
tol=1e-16;

Aligned = X;
G1 = 1e16; G0 = 1e15;

numIter = 0;

while (G1-G0) > tol
G1 = G0;
for i=2:k;
            [D, Z] = procrustes(Aligned(:,:,1), Aligned(:,:,i))
Aligned(:,:,i) = Z;
end
G0 = getsumsq(Aligned);
           numIter=numIter+1;
end


function ssq = getsumsq(X)
ssq = 0;
k=size(X,3);
for i=1:k-1
    for j=i+1:k
        ssq = ssq + sum( sum( ( X(:,:,i)- X(:,:,j) ).^2 ));
    end
end

"Shyamal L" <lshyamalNOSPAM@yahoo.com> wrote in message <iadqfl$7q7$1@fred.mathworks.com>...
> ........
> while (G1-G0) > tol
> G1 = G0;
> for i=2:k;
> [D, Z] = procrustes(Aligned(:,:,1), Aligned(:,:,i))
> Aligned(:,:,i) = Z;
> end
> G0 = getsumsq(Aligned);
> numIter=numIter+1;
> end
> ........
- - - - - - - - -
  I do not see the point in your iteration process. You are aligning each of the matrices Aligned(:,:,i) for i = 2:k with the single matrix Aligned(:,:,1) to the best of 'procrustes' ability - that is, to the matrix that possesses the least square differences with Aligned(:,:,1) via a translation, a unitary transformation, and a scale factor change. You do not change Aligned(:,:,1). Therefore once performed, you can do no better on subsequent passes, so why are you repeating these with your while-loop? If it repeats more than once, it will be due purely to round off differences. The G0 value will remain essentially unchanged after the first trip through except for these rounding differences.

  (By the way, "(G1-G0) > tol" is ordinarily not the way to check for differences outside a tolerance. You need "abs(G1-G0) > tol". However, in the above case it may have saved you from an infinitely repeating while-loop.)

Roger Stafford