siavashk / pycpd

Pure Numpy Implementation of the Coherent Point Drift Algorithm
MIT License
513 stars 115 forks source link

Question about tolerance #36

Closed gattia closed 4 years ago

gattia commented 4 years ago

I've been playing with speeding up PyCPD and have been using Cython to write a few functions and have gone back to the original paper / matlab code to implement some of the other methods (low rank in particular).

Anyways, in that process, I've noticed that you don't always use the same error functions (tolerance) as they used in the matlab scripts.

For example, your deformable registration uses the change in sigma between two iterations where as the CPD original paper calculates an objective function that I believe is the negative log likelihood and uses a change in this to determine when the algorithm has converged.

Just curious if you remember if there was a reason you altered the tolerance/error? Im finding that using the original implementation that the algorithm is "converging" too early and therefore isn't actually getting the right transformation - I don't pass the tests. Im using the tests you setup as a means of determining if my changes are breaking anything.

Thanks!

siavashk commented 4 years ago

I've been playing with speeding up PyCPD and have been using Cython to write a few functions and have gone back to the original paper / matlab code to implement some of the other methods (low rank in particular).

This is a good idea and I am interested in seeing your results.

For example, your deformable registration uses the change in sigma between two iterations where as the CPD original paper calculates an objective function that I believe is the negative log likelihood and uses a change in this to determine when the algorithm has converged.

Yes, you are right. I found this to be empirically a better way of determining convergence. I do something similar in my GMM-FEM implementation too.

Im using the tests you setup as a means of determining if my changes are breaking anything.

I did not think too hard about convergence when I was writing the code. If you can figure out a mathematically rigorous way of determining convergence that will be great. Most of my test cases are bare minimum sanity checks, they might not be the best way of determining correctness.

gattia commented 4 years ago

Thanks for the reply and for clarifying/confirming. I just wanted to make sure I was interpreting correctly. And while the tests may not be mathematically rigorous, if the program isn't returning the right transformation (and it previously was) then there is clearly something wrong. And from that standpoint, they've served their purposes at least for my playing around). If I make any changes/additions to them I'll definitely let you know.

If you're interested in some of my other changes I actually pushed another repo (cycpd; https://github.com/gattia/cycpd) so that I/a colleague could use the new version. I found that re-writing the E-step in Cython reduces computation times for the rigid/affine to be ~1/3 of their time with PyCPD - I actually wrote/include two Cython versions of the E-step one that is just an optimization of the PyCPD implementation, and another that uses the same mathematical logic that Andriy (original) used - Andriy's implementation is a tiny amount faster - maybe not worth it for the lost clarity of what is being done.

As for the lowrank, as of right now that is entirely still in python (at the same cycpd repo). lowrank makes a HUGE difference in computation time. For a 5k point mesh the M-step in PyCPD takes ~11s on my laptop, the lowrank (with 100 eigenvectors/eigenvalues) takes <0.1s. So, when you use lowrank, the E-step becomes the bottleneck for all 3 methods.

The last thing to implement from the original paper is the FGT which is used on the E-step and seems like it would considerably speed up all three algorithms further.

siavashk commented 4 years ago

This is all great. For the record, I have invited you to the repository. Feel free to add your low rank approximation to the development branch if you have time. Just make sure you only use numpy since pycpd is intended for educational use.

gattia commented 4 years ago

Excellent! Thanks for the invite. When I get some time I'll incorporate the low rank to PyCPD as well.

Just saw that you pushed the new changes and updated to v2.0.0. Nice! I was mostly matching the format of PyCPD on CyCPD (so that people might better be able to learn/understand/follow), I'll probably make a few changes there to match the new naming conventions etc. If you're interested in contributing anything to CyCPD let me know.

Thanks again for working on and maintaining this! Its been of great utility to me, and I've learned lots via it.