Procrustes functional alignment bug

[shaangao]

Hi nltools developers,

Thank you for all your great work in developing this package -- I enjoyed using it so far!

One issue I encountered is the Procrustes functional alignment algorithm (in function align() when method == 'procrustes') not working as expected. I reviewed the source code in reference to the pattern of outputs I got, and a potential place to look into is a function call to scipy's orthogonal_procrustes() in the procrustes() function (lines 1513-1514 in stats.py):

# transform mtx2 to minimize disparity
R, s = orthogonal_procrustes(mtx1, mtx2)
mtx2 = np.dot(mtx2, R.T) * s

In nltools' implementation of procrustes(), data1/mtx1 is supposed to be the reference, and data2/mtx2 is the data to be transformed with the objective of minimizing its Frobenius norm with data1.

According to scipy's documentation for orthogonal_procrustes(), however, the first argument (corresponding to mtx1 in the function call in line 1513) is the matrix to be transformed, and the second argument is the reference matrix. Therefore, could you try updating line 1513 (and other related places if any) to

R, s = orthogonal_procrustes(mtx2, mtx1)

and check if this will make procrustes() work as expected?

Thanks much! Shan

[ejolly]

Hey @shaangao thanks for the kind words and for noticing this. We actually use the generalized procrustes which is the same as scipy.spatial.procrustes() see docs, i.e. de-means and normalizes matrices first. In fact our implementation is identical (same order of arguments) except we pad with 0s if the matrices don't have the same shape.

It's a little confusing because procrustes and orthogonal_procrustes have the order of their arguments reversed, but procrustes makes a call to orthogonal_procrustes and uses the transpose of the solution to apply the transform.

This works because using the transpose of R is equivalent to swapping the order of arguments in orthogonal_procrustes():

R, s = orthogonal_procrustes(m1, m2)
R2, s = orthogonal_procrustes(m2, m1)

np.allclose(m2 @ R2, m2 @ R.T) # True

[ejolly]

Related to #417

[ejolly]

I'm going to close this for now @shaangao , but feel free to reopen if you have a related issue!