Re: [eigen] Progress and plans

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

2010/2/18 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> 2010/2/18 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
>> Hello,
>>
>> I thought it might be useful to give you a general update so that you can
>> take it into account when making plans during the meeting. I hope that
>> you'll be able to report some of your conclusions back to the list.
>>
>> The MatrixFunctions module is almost ready for release (more details further
>> down). I plan to continue working on Eigen, but the course I'm teaching at
>> the moment is proving a lot of work so I have little spare time. Thus I'm
>> glad I decided not to come to the meeting. The students are having their
>> Easter break starting 22 March; I hope that will give me a bit of air.
>> Initially, my plan was to work on precision-oriented testing and I thought a
>> bit about it (see below), but at the moment I'm not sure whether I'll have
>> enough time to have something in place in time for 3.0. An alternative would
>> be for me to improve the tutorial, which has the advantage that even a
>> little time yields results.
>>
>> I welcome suggestions on what areas it would be most useful for me to
>> contribute to.
>>
>>
>> Status of MatrixFunctions module
>> --------------------------------
>>
>> The module contains implementations for good algorithms for computing the
>> matrix exponential and (hyperbolic) sine and cosine, and entire functions in
>> general. The main function missing is the logarithm.
>>
>> I think the module is usable at the moment. However, there are some problems
>> if the matrices have clustered eigenvalues (this is the most difficult
>> case). One issue is that ComplexSchur often reaches the maximum number of
>> iterations. I don't know whether that's an inherent weakness in the algoritm
>> or only in the implementation - it may be as simple as increasing the
>> maximum number of iterations.
>
> It'd be interesting to know if increasing the number of iterations helps.
>
> I've been thinking about the following idea for this kind of iterative
> processes for diagonalization/Schur/SVD: we could add a termination
> condition of the form "the last 10 iterations haven't increased
> precision in a significant way". In addition, of course, to the
> condition "the matrix is within machine precision in the wanted form".
>
> At least I'll be experimenting it in SVD.
>
>> There seems to be another issue in
>> matrix_function_3, which I haven't investigated yet. My guess is that the
>> test matrices are simple too difficult to evaluate the matrix functions at
>> the desired precision, which is dummy_precision() for floats = 1e-5 = 200 *
>> machine epsilon. The test with double seems to go okay and dummy_precision()
>> for doubles = 1e-12 = 10000 * machine epsilon; any reason why we're so much
>> stricter with floats?
>
> No particular reason. But keep in mind that dummy_precision is just
> that, a dummy precision value for use when you don't want to think
> hard about the precision that you want. Also keep in mind that when,
> in unit tests, you use e.g. VERIFY_IS_APPROX, it uses another, yet far
> coarser precision level, given by test_precision() in test/main.h.
>
> In your case, it probably makes sense to use a custom precision level,
> doing something like
>
> VERIFY(result.isApprox(expected, my_precision));
>
>>
>> As far as I am concerned, the API is quite okay. The module provides
>> functions like
>>
>>   template <typename Derived>
>>   MatrixExponentialReturnValue<Derived>
>>   ei_matrix_exponential(const MatrixBase<Derived> &M)
>>
>> The return type is a derived class of ReturnByValue. Issues I'd like to
>> discuss:
>> * is the ei_ prefix necessary, given that the function is in the Eigen
>>  namespace? The wiki has the rule "global functions start with ei_" .
>
> Hm, we should update that: it applies mostly to function names that
> would be too polluting without the ei_ prefix. For exemple ei_sin and
> ei_exp. It can also serve to hint that a function is for internal use
> only. In your case, don't put a ei_ prefix. And since elsewhere we're
> abbreviating exponential as exp, you can do that too. So, I see two
> possible names for ei_matrix_exponential:
> 1) matrixExp()
> 2) ei_exp(). After all, we're defining it for numbers, so why not for
> matrices. Unfortunately, without ei_ it would conflict with std::exp
> for users who do "using namespace".
>
>> * should it be a member function instead? Is this possible while
>>  maintaining the separation between stable and unstable modules?
>
> That is perfectly possible. The prototype would have to be in
> MatrixBase.h, so in Core, but that's not a problem. Just because the
> prototype is there doesn't mean it's "stable": users still can't use
> it without explicitly including your module.
>
> So:
>  matrix.exp()
>
> That's certainly the most elegant possible syntax!
>
>> * related, what needs to be done to get the module out of unstable/ ?
>
> You're doing pretty much exactly what you have to :)
> 1) say that you want it out of unstable/ ---> check
> 2) make sure the API is ready for prime time ---> in progress
> 3) tests, documentation ----> check
> 4) commit to maintaining that stuff in the future ? ;)

The implicit thing here is that, in the case of your module, I have no
concern about its potential usefulness. Several people have requested
that already, and it's a reasonable feature.

Of course, that's only my opinion, it would be good to hear from the others ;)

Benoit


>
> Benoit
>
>>
>>
>> Precision-oriented tests
>> ------------------------
>>
>> I thought a bit about it, nothing profound, but perhaps you're interested.
>> The ToDo list on the wiki seems to indicate three possibilities:
>>
>> 1. We write an BLAS/LAPACK interface with Fortran (or C?) bindings so that
>> we can run the LAPACK test suite directly.
>>
>> 2. We port the LAPACK test suite to call the Eigen routines in C++.
>>
>> 3. We write our own test suite.
>>
>> My preference goes to option 2. It's difficult to find good test matrices,
>> so leveraging LAPACK's work is good. Furthermore, it's a good selling point
>> to be able to say that we pass the LAPACK test suite. My guess is that
>> option 1 is more work than 2, especially to make it portable. Fortran/C
>> bindings may be useful in themselves, but we will be up against a crowded
>> field.
>>
>> The LAPACK test suite is described in LAPACK Working Note 41, available from
>> http://www.netlib.org/lapack/lawns/lawn41.ps
>>
>>
>> Have a good meeting!
>> Jitse
>>
>>
>>
>