add procrustes?

[argriffing]

The problem is to find a transformation of n p-dimensional points X_i, maybe allowing translation/rotation/scaling, so that the transformed points closely match up with n corresponding target points Y_i with mismatch penalized by maybe sum of squares of errors. I've written code for this at one time, downstream projects e.g. skbio.math.stats.spatial.procrustes also have implementations, matlab has it, web forums have questions about how to do this in numpy/scipy, and it would seem to be in scope to be added to scipy.
http://en.wikipedia.org/wiki/Orthogonal_Procrustes_problem

[jacobcvt12]

I think this would be nice to add in to scipy. I haven't contributed to any scipy enhancements before, but I'd like to try this one/help out with it, if the decision is made to add in this feature.

[ElDeveloper]

If there's interest to add this to scipy, we would be willing to contribute the implementation in skbio.math.stats.spatia.procrustes.

cc @gregcaporaso, @jairideout, @ebolyen, @antgonza, @rob-knight, @wdwvt1

[antgonza]

Sounds good to me, if there is nothing there to do that. 👍

[rgommers]

Looks good to me. I see skbio puts it in spatial, but it looks like in scipy it would belong in linalg right?

Could you bring it up on the mailing list? That's where decisions on including new functionality are made.

[ElDeveloper]

Sounds good I am waiting for the e-mail confirmation from the scipy-dev list, once I get that I'll reference this issue there.

[sturlamolden]

I think it would be better to return the orthogonal transformation matrix omega. In case of a solver for the orthogonal Procrustes problem in scipy.linalg it should be the only result that is returned.

[sturlamolden]

Procrustes analysis is regarded a spatial statistics problem. But solving the orthogonal Procrustes problem is a linear algebra problem. I think a low level solver for the orthogonal Procrustes problem would fit best in scipy.linalg, but a high level Procrustes analysis (with full statistical tests) could go in scipy.spatial.

[sturlamolden]

Also, don't use np.transpose(a) in the code. Using a.T is preferable, and np.dot is smart enough to handle C and Fortran arrays correctly with transpose flags in BLAS.

In SciPy, code that uses scipy.linalg.svd should preferably use Fortran arrays consistently to avoid f2py overhead. So the procrustes solver would typically start with calling np.asfortranarray on its inputs. This also allows the input to be "array-like".

[argriffing]

@sturlamolden Maybe a scipy.linalg.orthogonal_procrustes(A, B) function which uses the Wikipedia notation A, B, and R, and which is implemented similarly to the two lines

u, s, vh = np.linalg.svd(np.dot(np.transpose(mtx1), mtx2))
q = np.dot(np.transpose(vh), np.transpose(u))

in the scikit-bio helper function?

[sturlamolden]

scipy.linalg.procrustes would be a shorter name.

I would also suggest thet we use scipy.linalg.svd instead of numpy.linalg.svd in SciPy.

The main thing would be that a function in scipy.linalg.* should solve the linear algebra problem, and nothing else. Cf. that scipy.linalg.lstsq solves the OLS problem, but does not actually do linear regression analysis.

scipy.spatial.procrustes could do full Procrustes analysis and call scipy.linalg.procrustes as a utility function.

[argriffing]

Added a PR #3809 for the orthogonal part that is purely linear algebra, without translation or scaling. This follows the API suggested by @sturlamolden, but with the longer function name.

[argriffing]

As I was adding that PR, I noticed small inefficiencies in the scikit-bio procrustes code which I thought I'd write down while I still remember. No extra dimension is needed for reflection. If you remove the extra-dimension-related code you will find that the 2d example in the docstring, which requires reflection, will still work. Also, code like np.sqrt(np.trace(np.dot(result, np.transpose(result)))) can be replaced by norm(result).

[ElDeveloper]

@sturlamolden, @argriffing should I wait for #3809 to be merged and then work on porting over the procrustes code with the improvements that @argriffing has pointed out thus far? Otherwise what would be the best way to do this?

[argriffing]

I tried setting up scikit-bio for development but its tests keep failing and I'm not patient enough to finish fixing it, so I probably won't make a scikit-bio PR for those improvements soon. Currently I think I have an incompatible version of "mpl_toolkits."

[gregcaporaso]

@argriffing, what issues were you having with scikit-bio install? pip install works for us on linux and mac (though we don't currently support windows, so we expect issues there).

[argriffing]

what issues were you having with scikit-bio install?

When I run nosetests --with-doctest skbio from my home directory as suggested on http://scikit-bio.org/ I get the error

[...]/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 1040, in cla
        self.zaxis._set_scale('linear')
    AttributeError: 'ZAxis' object has no attribute '_set_scale'

I assume this is because I screwed up the installation or in trying to set up a development environment or because I have an obsolete version of some package. This is on Ubuntu with Python 2.7. If I find myself motivated to try to fix it then if I have any questions I'll be sure to ask on a scikit-bio forum!

[argriffing]

Two skbio PRs are submitted!

[ElDeveloper]

Thanks!

On (Jul-22-14|18:15), argriffing wrote:

Two skbio PRs are submitted!

---

Reply to this email directly or view it on GitHub:
#3786 (comment)

[argriffing]

Just to give an update for this issue: the two skbio PRs are now merged, I think the scipy linalg procrustes PR should be in a reviewable state, and the hypothetical scipy spatial procrustes PR has not yet been created.

[ElDeveloper]

@argriffing, I just added one minor documentation comment to #3809! Thanks, I'll try to get a PR ready in the next few days.

[ev-br]

Linalg procrustes PR merged in e346aa5, can this be closed?

[argriffing]

@ev-br this isn't quite finished, because the merged orthogonal_procrustes just implements the core linear algebra common to all procrustes-related applications. It is not really useful by itself. I think there is interest in adding a spatial data analysis layer on top of it, so that if you pass two matched sets of points to a function, then the function will return enough information to let the caller easily transform a new point from the first space towards the second space, subject to some constraint. For example you could allow only transformations consisting of compositions of translation, rotation, reflection, and scaling, but not shear. On the other hand if it's bad form to have multiple github PRs associated with a single github issue then a new issue could be opened and this one could be closed.

[ev-br]

Ok, let's keep this open then.