Re: [eigen] Recursive Matrices in Eigen

[Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx>]

Hi!

On 29/05/2019 10.21, Darius Rueckert wrote:
> Hello everyone,
>
> for the past months we have been working on making recursive matrices 
> work in eigen. Here is an example:
>
> SparseMatrix<Matrix<double,4,3>> A;
> Matrix<Matrix<double,4,1>,-1,1> x,b;
> x = A * b;

We've had some more or less stable support for
  [Sparse]Matrix<Array<...> >
for a while: 
https://bitbucket.org/eigen/eigen/src/3cbfc2d7/Eigen/src/Core/NumTraits.h#lines-233

Of course not quite the same, and probably adding sub-matrix support is 
slightly more complicated than sub-arrays (arrays behave much more 
similar to scalars than matrices).

> The obvious advantage is that we get support for block matrices on all 
> eigen types for "free" ( DiagonalMatrix,TriangularView,LDLT,...). After 
> some testing, the performance appeared to be "better than expected" and 
> we have decided to write a paper on how to apply recursive matrices on 
> non-linear least squares problems. The paper is accepted and will be 
> published at HPG in July. Here is the current version:
>
> https://cloud9.cs.fau.de/index.php/s/UINyrgAUtaAiWvR

That's nice! (Only glanced through your paper so far)

> Unfortunately there are two limitations which I was not able to fix.
>
> 1. All inner matrix types require an additional wrapper class, which I 
> have called "MatrixScalar". The matrix A from above is then:
>
> SparseMatrix<MatrixScalar<Matrix<double,4,3>>> A;
>
> 2. Some eigen kernels do not work directly. For example, some operations 
> use the function scaleAndAddTo(..., Scalar(1)), which does not work if 
> "Scalar" is a recursive matrix type. To fix this I had to specialize a 
> lot of eigen kernels, which creates redundancy with only minor changes.

I guess replacing `Scalar` by `NumTraits<Scalar>::Literal` here could 
solve the issue. `Literal` needs to be a typedef which can be implicitly 
multiplied by `Scalar`. This could also make code for some custom 
Scalars much more efficient (e.g., multiplying an AutoDiff type by a 
double, does not really require to first convert the double to the 
AutoDiff type). The main reason we are converting numbers to `Scalar` is 
to be compatible with custom scalars which do not allow implicit conversion.

> So here are my questions to you:
>
> What do you think about this idea?
>
> Do you think a deep integration into eigen would be possible? (Without 
> recursive specializations and without "MatrixScalar" class)

I generally like the idea. Of course the `MatrixScalar` intermediate 
should be avoided (I did not look into details yet why you need this).
Ideally, this should be possible with only adding a few Traits 
specializations (and probably a lot of cleanup like using `Literal` as 
suggested above).

I can't guarantee having much time soon for deeper discussions, but feel 
free to ask technical questions (perhaps make a pull-request for that, 
once it gets too technical).

> The code can be found here:
> https://github.com/darglein/EigenRecursive
>
> Example for sparse block-matrices:
> https://github.com/darglein/EigenRecursive/blob/master/samples/sparseOperators/sparseOperators_main.cpp 
>
>
> With best regards,
> Darius

--
 Dr.-Ing. Christoph Hertzberg

 Besuchsadresse der Nebengeschäftsstelle:
 DFKI GmbH
 Robotics Innovation Center
 Robert-Hooke-Straße 5
 28359 Bremen, Germany

 Postadresse der Hauptgeschäftsstelle Standort Bremen:
 DFKI GmbH
 Robotics Innovation Center
 Robert-Hooke-Straße 1
 28359 Bremen, Germany

 Tel.:     +49 421 178 45-4021
 Zentrale: +49 421 178 45-0
 E-Mail:   christoph.hertzberg@xxxxxxx

 Weitere Informationen: http://www.dfki.de/robotik
  -------------------------------------------------------------
  Deutsches Forschungszentrum für Künstliche Intelligenz GmbH
  Trippstadter Strasse 122, D-67663 Kaiserslautern, Germany

  Geschäftsführung:
  Prof. Dr. Jana Koehler (Vorsitzende)
  Dr. Walter Olthoff

  Vorsitzender des Aufsichtsrats:
  Prof. Dr. h.c. Hans A. Aukes
  Amtsgericht Kaiserslautern, HRB 2313
  -------------------------------------------------------------