yade-dev team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Bruno Chareyre said:     (by the date of Fri, 04 Jun 2010 16:28:35 +0200)

  this->call_count++;
  if(call_count % 10 != 0)
    return; // check only every 10th time
  static const Real MAX_ERROR=200*std::numeric_limit<Real>::epsilon();
  static const Real MIN_ERROR=3*std::numeric_limit<Real>::epsilon();
  Real sq_length(a*a+b*b+c*c+d*d)
  if(std::abs(sq_length - 1.0) > MIN_ERROR )
  {
    Real length = std::sqrt(length);
    if( std::abs(length - 1.0) < MAX_ERROR )
    {
      a/=length;
      b/=length;
      c/=length;
      d/=length;
    }
    else
      throw or exit(1) or something();
  }
  return; // OK

disclaimer: I typed that code without looking eigen's implementation
of quaternions.


The numerical error in quaternion becomes significant after roughly
10 to 20 multiplications of a given quaternion.


I say: check & normalize every 10th non-const operation on
quaternion. Do this in the library, in the quaternion class member
functions, not in user code.


> There is something I don't get. You (Janek) say we don't need to 
> normalize all the time, it is enough to check if length=1, but it is 
> almost the same. In these two lines the biggest cost is to compute the 
> norm :

you can check just by calculating squared norm, which is faster. The
most expensive operation is square root. Later to divide you need to
perform that square root operation.


> Real norm=q.norm(); if (norm>1+epsilon || norm<1-epsilon){
> q/=norm;//or whatever the function}

better to do it in a way that I've shown above: defer calculating
square root until later stage.

> (Janek) say we don't need to 
> normalize all the time

yes. exactly. Our quaternion is not unitary only if it diverges
from 1.0 more than THREE times epsilon(). That THREE comes from three
additions in this code: a*a+b*b+c*c+d*d. The epsilon() here is
potentially added here three times. And it is possible that a certain
combination of numbers leads to situation in which normalizing
produces numbers that are "just below", and you could normalize the
quaternion infinitely and never reach 1.0 (always 3*epsilon()
difference), because of numerical precision.

I didn't make that 4*3=12 times epsilon(), because we want to be on
the safe side. Better to normalize than not, just do not normalize if
that doesn't give any improvement in quality.

But OTOH if we reach error bigger than 200*epsilon() then it is
impossible to recover, and calculations must be interrupted.

There is another culprit: if we want to go faster by normalizing
every 10th call, then we need to add another variable to quaternion
class! It would be a counter of calls, and would count from 0 to 10.
The counters are independent for each quaternion. That certainly has
to be a private variable in class implementation. I'm not liking
it, but I don't see any other solution for someone who is concerned
with speed. We could benchmark to see if that "every 10th call" is
any significant change in the calculation speed.

I don't foresee myself implementing that in eigen quaternions anytime
soon. I hope that Vaclav or you will jump in ;) I will certainly
review your commit concerning this.

best regards
Janek Kozicki

> It is at least good to hear from you guys that unit quaternions still 
> give unit quaternions (put aside numerical errors) in eigen as in wm3.
> 
> Bruno
> 
> 
> Janek Kozicki a écrit :
> > Hi,
> >
> > you shouldn't be running wild through the code adding
> > quaternion.normalize() every second line in the code.
> >
> > I've went through this stage already, and let me explain what is the
> > problem and how to proceed.
> >
> >
> > Quaternions:
> >
> > Quaternions, are just like complex numbers, except that instead of
> > single imaginary "i", they have three imaginary numbers: "i", "j", "k".
> >
> > We are using them to perform rotations in 3D space. Namely, we are
> > using them as a representation of an SO(3) rotation group:
> > http://en.wikipedia.org/wiki/SO(3). Fortunately multiplying the
> > quaternions is just like composition of two rotations.
> > (BTW: orthogonal rotation matrices 3x3 are also a representation of SO(3) ). 
> >
> > Quaternions, are quite like complex numbers, and can be of ANY
> > magnitude (length) that you could imagine. You can calculate sin(q),
> > ln(q), e^q, solve polynomial equations with them, and other stuff.
> > Which we are not interested in.
> >
> >
> > The problem:
> >
> > The problem is that any given quaternion belongs to the SE(3)
> > rotation group only, and only, if its length equals ONE.
> > If any quaternion suddenly doesn't have length=1, it does not belong
> > any more to SE(3).
> >
> > Analytically multiplying quaternions with length=1 will always
> > produce a quaternion with length=1.
> >
> > Numerically multiplying two quaternions requires 12 additions and 16
> > multiplications. There is a numerical error of order 1e-17 which can
> > build up, making the resultant quaternion<double> to be not unitary
> > any more.
> >
> >
> > The way to solve this:
> >
> > Analytically + theoretically if a quaternion's length is not 1, then
> > this is for us a critical error since it is NOT a rotation any more.
> >
> > Numerically + practically there is a numerical error introduced by
> > the computer. It is of magnitude std::numeric_limit<TYPE>::epsilon()
> >
> > http://www.cplusplus.com/reference/std/limits/numeric_limits/
> > http://msdn.microsoft.com/en-us/library/6x7575x3%28VS.71%29.aspx
> >
> > (where used TYPE is usually float, double or long double, for double
> > IIRC it was 1e-17).
> >
> > We have to cope with this error, and due to this error it is
> > impossible to have a quaternion's length always equal to 1.
> >
> >
> > How this epsilon() propagates through calculations (12 additions, 16
> > multiplications) depends exactly on the way in which it was
> > implemented. The worst case would be: a single epsilon() becomes
> > 12*16*epsilon() larger after single multiplication. I am sure that
> > this is never the case. But you get the picture. We can call that
> > 12*16 a MAX_ERROR and use 200 instead of 12*16=192.
> >
> >
> > The problem again:
> >
> > So, we want to use UNITARY quaternions, but the library (wm3 and
> > eigen) provide quaternions ONLY. Not the unitary quaternions.
> >
> > Our problem is inherent to the mathematics of SO(3) and computer's
> > numerical errors. The end user must never worry about it.
> >
> >
> > The solution:
> >
> > We need a library that provides UNITARY quaternions.
> >
> > If we don't have such library, then we need to normalize quaternions
> > every second line in the code. And suddenly YADE's code base is
> > bigger by 5000 lines.
> >
> > The only correct way to cope with this, is to extend existing
> > quaternion library with a simple function that checks if quaternion
> > is UNITARY.
> >
> > if quaternion's length differs from 1, then:
> >
> > - if it differs less that MAX_ERROR*epsilon(). Then the error is not
> >   critical. We have to normalize that quaternion (divide it by its
> >   length) and silently move on.
> >
> > - if the difference from 1 is bigger than MAX_ERROR*epsilon()
> >   then YADE has to report a critical error, throw exception and instantly
> >   terminate all calculations. Since whatever calculations there were,
> >   they ARE wrong.
> >
> > This simple function has to be called at the end of each non-const
> > member function in quaternion class.
> >
> >
> >
> > Let me repeat the solution:
> >
> > 1) Do NOT run wild through the code normalizing quaternions everywhere.
> >
> > 2) Do check if quaternion is unitary at the end of each non-const
> > member function of quaternion class.
> >
> > 3) quaternions are in eigen lib and we can not modify it? Then we
> > have to either:
> >    - send them a patch, witch enables the library to provide unitary
> >      quaternions, or
> >    - derive from their class and extend it.
> >
> >
> > I hope that it is clear now.
> >
> > This error is NOT of type: OMG all our calculations were wrong!
> >
> >
> > The error build-up exists, usually is negligible. It shows itself only if:
> >
> > - same quaternion is used a lot in a single loop. lading to lots
> >   of multiplications in single iteration. This is not the case in
> >   DEM. Such multiplication in DEM with spheres is done just once per
> >   iteration per quaternion.
> >
> > - a quaternion is mistakenly initialized to be non-unitary. Which CAN
> >   be the case quite often, if you don't read what you are writing.
> >
> >
> > I stumbled on this error with lattice model, because rotating a
> > single rod required about 10 to 20 quaternion multiplications per
> > iteration per quaternion. And I thought that it should be fixed in miniwm3, but:
> > - I was too lazy to modify miniwm3
> > - I was scared of computation overhead that this is going to introduce
> > - I saw that in fact this problem manifests itself rarely.
> >
> >
> > In fact that function instead of doing unitarity check on every call,
> > could count the number of times that it was called, and perform
> > normalization only every 10th call. That solution doesn't seem to be
> > "perfect", but I tend to like it.
> >
> > Maybe it is better to do all those checks ALWAYS (each time) in DEBUG mode.
> >
> > And do them every 10th time in non debug mode. You cannot skip them
> > totally in non-debug mode, because in non-debug mode the numerical
> > buildup error occurs too!
> >
> >
> > well, I hope that now you understand what's going on here, and aren't
> > scared anymore
> >   
> 
> 
> -- 
> _______________
> Bruno Chareyre
> Associate Professor
> Grenoble INP
> Lab. 3SR
> BP 53 - 38041, Grenoble cedex 9 - France
> Tél : 33 4 56 52 86 21
> Fax : 33 4 76 82 70 43
> ________________
> 
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help   : https://help.launchpad.net/ListHelp
> 


-- 
Janek Kozicki                               http://janek.kozicki.pl/  |