Re: [eigen] Recursive Matrices in Eigen

[ Thread Index | Date Index | More Archives ]


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:

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:

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:

Example for sparse block-matrices:

With best regards,

 Dr.-Ing. Christoph Hertzberg

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

 Postadresse der Hauptgeschäftsstelle Standort Bremen:
 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:
  Deutsches Forschungszentrum für Künstliche Intelligenz GmbH
  Trippstadter Strasse 122, D-67663 Kaiserslautern, Germany

  Prof. Dr. Jana Koehler (Vorsitzende)
  Dr. Walter Olthoff

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

Mail converted by MHonArc 2.6.19+