Re: [eigen] Arbitrary precision arithmetic support in Eigen

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

2010/7/15 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> On Thu, Jul 15, 2010 at 8:37 AM, Pavel Holoborodko
> <pavel@xxxxxxxxxxxxxxx> wrote:
>> Hi Gael!
>>
>> I am very thankful for your warm acceptance of MPFR C++.
>> Thank you for your support!
>>
>> 1.
>> I've added 3 new functions to MPFR C++ library:
>> mpfr::machine_epsilon(prec) // the smallest eps such that 1+eps != 1 for the
>> precision prec
>> mpfr::mpreal_min(prec)         // the largest positive number presentable in
>> precision prec and current exponent range
>> mpfr::mpreal_max(prec)        // the smallest positive number presentable in
>> precision prec and current exponent range
>>
>> They do not use pow, only simple and fast operations.
>> Also I've changed mpreal_eigen.h to embed updated MPFR C++.
>
> excellent! seems to work well. I'm currently only testing high-level
> operation such LU, LLT, and symmetric eigenvalues, that indirectly
> already covers a large part of Eigen !!
>
>
> To Eigen people: maybe this is an opportunity to add new accuracy
> oriented unit tests that would compares the results obtained with
> float or double to the ones obtained with a much higher accuracy?

I'm not sure that this is the right approach to precision tests:

First, we can test precision just as well without using higher
precision numbers, by computing something whose result is well know
from the input data. For example, if we want to test the precision of
a decomposition, we can simply reconstruct the matrix from the
decomposition and compare with the original matrix.

Second, using bignums will be slow, so in practice we'd have either
very slow tests or we'd be giving up testing large sizes, which is
important in precision tests.

Third, certain kinds of floating-point precision issues might persist
even with much higher precision, so just because a result agrees with
another computed with some higher precision, doesn't mean it's
correct. Of course if you keep increasing the precision, eventually
you catch all FP issues, but it's not clear beforehand how much you
need to increase it for a particular test. For example, if the issue
is that we do something stupid like

     (1 + x) - 1

and as a result we get 0 instead of the expected x when x = 1e-20 with
T=double, then for sure using quad precision will catch that, but it
won't catch the same issue with 1e-50, etc... just pointing out how
this kind of approach can fail to find bugs.

Benoit


>
>>
>> MPFR C++:
>> http://www.holoborodko.com/pavel/wp-content/plugins/download-monitor/download.php?id=4
>>
>> mpreal_eigen.h:
>> http://www.holoborodko.com/pavel/wp-content/plugins/download-monitor/download.php?id=9
>>
>> I would appreciate any tests, or further suggestions.
>>
>> 2.
>> I am very glad you fixed bug in ei_aligned_delete().
>> I believe that was the bug users complain about in KDE forum.
>> (crash on memory deallocation for dynamic matrices using mpreal)
>> http://forum.kde.org/viewtopic.php?f=74&t=43643.
>
> indeed, that was exactly it.
>
>> 3.
>> On where to put such a header file.
>> I agree with you that it would be convenient for users to include
>> MpfrSupport.h into Eigen itself.
>> Please do it.
>
> ok, done :) To be more precise I've currently added it as a
> MPRealSupport module which currently lives in the unsupported/ folder
> to not distract us too much for Eigen3. I've also added a
> FindMPFR.cmake file and to ease testing I've included a copy of
> mpfrc++ (the 2 files) in our test folder.
>
> Since mpfrc++ it is quite lightweight and not packaged in common
> distro, I'm even tempted to put these two files in the MPRealSupport
> module. We would update them every new release of mpfrc++. This way we
> can guarantee that the trio  MPFRC++/MPRealSupport/Eigen always work
> fine, and people would only have to install MPFR and #include
> <Eigen/MPRealSupport> and basta! Of course if you're OK for that we'll
> make it crystal clear that's come from MPFRC++ with link to your
> website everywhere it is appropriate (website, documentation, etc.).
>
> gael
>
>> One issue is maintenance of MpfrSupport.h upon changes in MPFR C++.
>> I think we can manage to do that properly and timely.
>>
>> Thank you! And please keep me updated.
>>
>> All the best,
>> Pavel Holoborodko
>> http://holoborodko.com/pavel/
>>
>> On Thu, Jul 15, 2010 at 5:27 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
>> wrote:
>>>
>>> while trying to write proper unit tests, I hit the following issue:
>>>
>>> mpfr::pow(2,-mpfr::mpreal::get_default_prec())
>>>
>>> which is used to compute the epsilon takes forever... Is it "normal" ??
>>>
>>> The following runs fine though:
>>>
>>> mpfr::mpreal x = 1;
>>> for(int i=0; i<mpfr::mpreal::get_default_prec(); ++i)
>>>  x = x/2;
>>>
>>>
>>> Also is there an easy way to get the highest and lowest representable
>>> values ? They are needed for some operations in Eigen3.
>>>
>>> cheers,
>>>
>>> gael.
>>>
>>> On Wed, Jul 14, 2010 at 9:24 PM, Gael Guennebaud
>>> <gael.guennebaud@xxxxxxxxx> wrote:
>>> > ok, I could not resit to adapt it for Eigen3. Actually this is very
>>> > useful to stress Eigen with scalar types having non trivial ctor/dtor
>>> > (i.e., allocating/freeing memory), and bingo (!) there was a bug in
>>> > ei_aligned_detete() for null pointers ;)
>>> >
>>> > gael
>>> >
>>> >
>>> > On Wed, Jul 14, 2010 at 4:06 PM, Gael Guennebaud
>>> > <gael.guennebaud@xxxxxxxxx> wrote:
>>> >> Hi Pavel,
>>> >>
>>> >> thank you for your email. This is excellent!
>>> >>
>>> >> Actually, the possibility to use arbitrary precision scalar types is
>>> >> one of the main argument I often use to sell Eigen. So having an
>>> >> officially supported module for your lib with extensive testing is
>>> >> something I wanted to do for a while, but you know there are always so
>>> >> many things to do.
>>> >>
>>> >> So I'm really glad to see you added support for Eigen ;) I will add a
>>> >> word on the website.
>>> >>
>>> >> What to do next?
>>> >>
>>> >> You said there is a lack of extensive testing. Since all our efforts
>>> >> are now focused on the next Eigen version, Eigen3 which should be
>>> >> release within a few months I'd rather spend time testing MPFR support
>>> >> on Eigen3 than Eigen2. However, the way to add support for custom
>>> >> scalar types slightly change. Fortunately, based on your current
>>> >> support header file this should not be difficult. I can have a look at
>>> >> it myself.
>>> >>
>>> >> One leading question is where to put such an header file? Within Eigen
>>> >> as a (MpfrSupport module) or as a standalone file as it is currently
>>> >> the case? The main advantage I see in having it in Eigen itself is
>>> >> that should simplify to had full and continuous testing What do you
>>> >> think?
>>> >>
>>> >> Finally, once MPFR C++ will use expression template, adding a good
>>> >> support for Eigen-MPFR C++ that leverages ET at both the matrix and
>>> >> scalar type level will be more tricky because Eigen will have to know
>>> >> about the returned expression types of the mpreal operations. We
>>> >> currently do not have an easy solution for that but this is something
>>> >> we have to workout anyway for many other use cases (e.g., to nest
>>> >> Eigen's Matrix/Array types, for automatic differentiation, etc.)
>>> >>
>>> >> cheers,
>>> >>
>>> >> gael
>>> >>
>>> >>
>>> >>
>>> >>
>>> >>
>>> >> On Wed, Jul 14, 2010 at 7:01 AM, Pavel Holoborodko
>>> >> <pavel@xxxxxxxxxxxxxxx> wrote:
>>> >>> Hello Everyone!
>>> >>>
>>> >>> First of all I would like to thank all the developers for the this
>>> >>> beautifully designed and implemented library!
>>> >>> I've learned a lot looking over it.
>>> >>>
>>> >>> I'm developer of C++ interface for the arbitrary precision floating
>>> >>> point
>>> >>> arithmetic library MPFR: http://www.holoborodko.com/pavel/?page_id=12.
>>> >>>
>>> >>> During last year I've received couple of questions from the users on
>>> >>> how to
>>> >>> use my wrapper together with Eigen. There was discussion in KDE forum
>>> >>> on the
>>> >>> topic but it seems that not all problems were solved:
>>> >>> http://forum.kde.org/viewtopic.php?f=74&t=43643.
>>> >>>
>>> >>> Recently the most complete (hopefully:)) bridge for embedding mpreal
>>> >>> type
>>> >>> into Eigen was created.
>>> >>> Please check home page under "Software using MPFR C++" for header file
>>> >>> mpreal_eigen.h:
>>> >>> http://www.holoborodko.com/pavel/?page_id=12.
>>> >>>
>>> >>> So, now user can equip Eigen with arbitrary precision arithmetic by:
>>> >>> 1. Installing  MPFR and GMP libraries.
>>> >>> 2. Using MPFR C++ interface (mpreal.h, mpreal.cpp) along with header
>>> >>> file
>>> >>> mpreal_eigen.h
>>> >>> 3. Using mpreal type (provided by MPFR C++) instead of built-in type
>>> >>> double
>>> >>> or float.
>>> >>>
>>> >>> Precision (number of bits used for significant) can be setup globally
>>> >>> by
>>> >>> mpreal::set_default_prec().
>>> >>> For instance double has only 53 bits precision with machine epsilon
>>> >>> 1e-16
>>> >>> (maximum 16 correct digits).
>>> >>> Using mpreal user can do calculations with any precision, e.g. 1024
>>> >>> bits
>>> >>> with machine epsilon 1e-309 (309 correct digits).
>>> >>>
>>> >>> I would appreciate if someone who intensively use/develop Eigen would
>>> >>> test
>>> >>> how MPFR C++ works with Eigen.
>>> >>>
>>> >>> Any suggestions for improvement, test reports, etc.  are very very
>>> >>> welcome.
>>> >>>
>>> >>> Also I would appreciate if Eigen webpage or documentation will
>>> >>> reference
>>> >>> such possibility on using arbitrary precision arithmetic.
>>> >>> I think, this option could be very useful for the users.
>>> >>>
>>> >>> Besides reference on my site would give me huge support and motivation
>>> >>> for
>>> >>> further improvement.
>>> >>>
>>> >>> All the best,
>>> >>>
>>> >>> Pavel Holoborodko.
>>> >>> http://holoborodko.com/pavel/
>>> >>>
>>> >>> P.S.
>>> >>> I'm planning to release new version of MPFR C++ soon based on
>>> >>> expression
>>> >>> templates inspired by Eigen.
>>> >>>
>>> >>>
>>> >>
>>> >
>>>
>>>
>>
>>
>
>
>