Re: [eigen] Geometry module - Quaternion fitting alternative

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/5/27 Hauke Heibel <hauke.heibel@xxxxxxxxxxxxxx>:
> On Tue, May 26, 2009 at 4:28 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> 2009/5/26 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>> wow that's a very great contribution !
>>>
>>> I did not read the code but I still have a few comments.
>>>
>>> On Tue, May 26, 2009 at 3:31 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>>>> 2009/5/26 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>>>>> This means that a prerequisite to vectorize that is that your matrices
>>>>> src and dst be RowMajor.
>>>>
>>>> Oh but I overlooked the fact that dst and src are the arguments passed
>>>> to the function. So we don't control their storage order, the user
>>>> does, and the default is col-major. We could say that you change the
>>>> implementation so that now the _rows_ are the points, so vectorization
>>>> can be done in the col-major case... but that would mean interleaved
>>>> storage, so the construction of the matrices in the first place would
>>>> be less efficient (less memory locality) so it's probably not worth
>>>> it.
>>>
>>> then I would say use the highest level Eigen expression as you can,
>>> and hope that in some situation the vectorization will be enabled. For
>>> instance there are still a few possibilities to improve the
>>> vectorization of partial reduction. Also note that the geometry module
>>> already requires the Array module.
>
> Well, here are right away a few important things. First, I Eigen'fied
> the covariance computation. The mean and the std deviation are now
> implemented purely based on rowwise(). The only thing which is
> currently done manually is the de-meaning.

OK. the next step is to add the "broadcasting" functionality to
PartialRedux so you can do rowwise()+=vector.

>
> Regarding the storage type (i.e. RowMajor vs. ColMajor) I would
> refrain from starting to implement any algorithm such that it makes
> any assumptions on it.

We agree with this but sometimes it's very hard to make code that will
be equally fast in both cases. Look at the PartialLU: it could equally
well be done with row operations or col operations, and of course
better speed is achieved if this matches the storage order. Currently
it is done with col operations and I'm not sure how to best make it
work in the row-major case without writing the code twice. Or is this
inevitable? Just transposing doesn't cut it as that would interchange
L and U.transpose() and they're not interchangeable (L has 1's on the
diagonal).

> In my programs I always default the storage
> type to RowMajor. I don't want to start a discussion here since I am
> convinced that the amortized cost in a typical program will be
> identical for both storage types and that the choice is merely a
> matter of taste or maybe interoperability in case you want to
> interface with other libraries.

Sure we agree with this. And to begin with, interoperability with
Eigen's own default which is, col-major, and therefore what a majority
of people use. So, given the choice, we take the path that goes faster
for col-major, but of course this isn't quite satisfactory and ideally
we should support both equally.

Since you're interested in Row-Major let me mention that by our unit
tests test mostly only col-major matrices so a very useful change
would be to change some matrix types to be row-major. For example if a
test tests MatrixXf,d,cf,cd, you could change the 'f' or the 'd' to
row-major. I say 'f' or 'd' more than 'cf' or 'cd' because one tricky
issue with storage order is vectorization and packet access, so it's
more interesting to test row-major on a type that is currently getting
vectorized.

>>> Also it would be very cool to describe the underlying algorithm in 2/3
>>> sentences so that the user know what to expect (complexity, etc.). If
>>> this is already the case just forgot what I said !
>
> I did that and added O-notation run-time approximations - it would be
> cool if someone could double check the values I have given. The most
> costly operation of the covariance computation (the umeyama?) is the
> matrix product of the demeaned input vectors - it's a product of an
> m-by-n times n-by-m matrix, where m is the dimension and n the number
> of input points.

Since m is typically small, that is not "much" more costly than the
rest, so this is a complicated algorithm to optimize: there's no
single big "critical part".

>
>>> I'm wondering whether the Geometric module is the right place for that
>>> ? More generally I'm wondering whether we should favor few but fat
>>> modules, or many but light ones ? For instance, if we want to go for
>>> the few but fat solution, then we could also have a single
>>> Decomposition module with LU, LLT, QR, SVD,, etc. !! Personally, I
>>> prefer the "many but light" solution. This is why I'm suggesting to
>>> put this algorithm in a new module, just like a Levenberg-Marquat
>>> solver should have its own module. Again, this is just a suggestion...
>>
>> Actually I 100% agree.
>>
>> Maybe the modules should be named after what they offer (like
>> "Umeyama" and "GeometricPolarDecomposition") rather than how they do
>> it (like "GeometryUsingSVD") so as to keep us freedom to reimplement
>> in the future.
>
> This is a huge topic - I say so because it is so tricky to decide
> about the membership. Here, I would still vote for the Geometry module
> because actually the algorithm is computing the geometrical
> relationship of two point sets. We could also put it in a module
> called 'Registration' because we could see the outcome of the
> algorithm as an alignment (registration) of data sets. Modules as
> 'Umeyama' and 'GeometricPolarDecomposition' sound so special to me
> that I am afraid of having like dozens of modules after a short period
> as soon as we follow that path. For Levenberg-Marquat it seems easier
> - possible something like 'Optimization'. I think in general it is
> more convenient for the user to find things when the modules are
> partitioned in a "What is the algorithm doing"-fashion as opposed to
> "How is the algorithm doing it"-way. Consider the case, when somebody
> pops up and implements the quaternion fitting. If that would end up in
> a different module than the Umeyama I would find it rather awkward
> considering that the algorithms are more or less trying to achieve the
> same.

I think we all agree on this but you enlightened me with more precises
examples than I had thought of. At this point I don't see a very
obvious way forward, maybe we can keep things as-is until we have an
idea...

>
> That said, for the initial commit I kept the algorithm in the Geometry
> module - until we (or you) have decided.

Sure.

Cheers,
Benoit