Re: [eigen] Progress and plans

[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 ? ;)

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
>
>
>