Re: [eigen] Allocation policy of eigen decompositions

On Fri, Mar 5, 2010 at 4:01 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

2010/3/5 Benoit Jacob <jacob..benoit.1@xxxxxxxxx>:

> 2010/3/5 Márton Danóczy <marton78@xxxxxxxxx>:
>> Dear Adolfo,
>>
>> and while you're at it - it would be also useful to have the
>> decompositions (or, at least in my case, HouseholderQR) work
>> optionally in-place. Internally, HouseholderQR copies its argument
>> into its member m_qr of the same type and performs the decomposition
>> on that, so it would be easy to expose this behaviour as an option, at
>> least for decompositions without pivoting. Benoit, you being the
>> master of the decompositions API, what do you think?

Thank you for putting that on the table because this is something I wanted to be discussed before 3.0 since it might potentially slightly affect the current API. Indeed, to make think more complicated, I'd like that the inplace API can also work on blocks... Yes, there are many use cases for blocks, and for them working out of place is a no sense!

In the meantime, if we forget about blocks, there is an (almost) already working solution. Since each decomposition has a ctor with only sizes to initialize the internals (maybe not all yet, but they should), then one can already do:

FullPivLU<MatrixType> lu(rows,cols);
MatrixType& m = lu.matrixLU(); // currently you need a const cast

for()
{
....
m = bla;
lu.compute(m); // the decomposition happens in-place :) (we should ensure that m is not copied into itself)
....
}

That works, but I agree that's not ideal because:
- you must declare the decomp first
- you cannot share the matrix with other decomp
- this does not address blocks

> Indeed that would be a great addition. Fortunately, here again, it
> seems to me that that can be done at any later date without impacting
> the existing API.
>
> Possible approach: add a computeInPlace() method, where you put 99% of
> the code of the existing compute() method, and rewrite compute() to
> just copy the matrix to a temporary and call computeInPlace().
>
> This is just adding API, doesn't require changing existing API. So
> again --- patches welcome, not my own priority :)

Wait! I said something stupid. All other methods in the decomposition
will still try to address the stored matrix member that never was
initialized. Bad!

A quick and dirty solution is to use: FullPivLU<MatrixXd&>

But that's not very good as it instantiates again a second time all
the code from FullPivLU<MatrixXd>.

The real solution is probably to store in addition to the member
matrix, a member pointer to a matrix.

All methods would use the pointer, not directly the member matrix.

The in-place decomposition method would set the pointer to point to
the matrix passed to it.

The non-in-place method would copy the matrix into the member matrix
and let the pointer point to that.

This sounds reasonale but does not address the cases of blocks. For blocks I have some ideas, I will come back to it later.

gael

Opinions?

Benoit

>
> Benoit
>
>
>>
>> Marton
>>
>> 2010/3/5 Adolfo Rodríguez Tsouroukdissian <dofo79@xxxxxxxxx>:
>>> Hello,
>>>
>>> I have a question on what the allocation policy of the eigen 3..x matrix
>>> decompositions will be, for this will strongly bias how I use eigen on
>>> different parts of my code.
>>>
>>> On the one hand, I know that the template parameters of the matrix to be
>>> decomposed (size, options, max size, ...) will be propagated to all the
>>> decomposition's members and temporaries [1], so if the input matrix is
>>> stack-allocated, so will be the internals of a decomposition. So far so
>>> good.
>>>
>>> On the other hand, I had the expectation that if heap-allocated matrices are
>>> used, all heap allocations of the decomposition algorithm would take place
>>> during the initialization phase, and none would be done in the computation
>>> phase (provided that the initialization and computation matrices have the
>>> same size in memory). This currently does not seem to be the case. For
>>> example,
>>>
>>> MatrixXd I(MatrixXd::Identity(16, 16));
>>> MatrixXd A(MatrixXd::Random(16, 16));   // Same size as I
>>>
>>> Eigen::FullPivLU<MatrixXd> lu(I);       // Allocates 5 temps
>>> lu.compute(A);                        // Allocates 2 extra temps
>>>
>>> Eigen::JacobiSVD<MatrixXd> svd(I);      // Allocates 4 temps
>>> svd.compute(A);                         // Allocates 1 extra temp
>>>
>>> I haven't checked where the extra allocations in compute() come from, but my
>>> question is, can and will they be prevented for 3.x? or is this an
>>> intentional design decision? (and if so, for what reasons).
>>>
>>> TIA,
>>>
>>> Adolfo
>>>
>>> [1]
>>> http://listengine.tuxfamily..org/lists.tuxfamily.org/eigen/2010/03/msg00095.html
>>>
>>>
>>> --
>>> Adolfo Rodríguez Tsouroukdissian, Ph. D.
>>>
>>> Robotics engineer
>>> PAL ROBOTICS S.L
>>> http://www.pal-robotics.com
>>> Tel. +34.93.414.53.47
>>> Fax.+34.93.209.11.09
>>> AVISO DE CONFIDENCIALIDAD: Este mensaje y sus documentos adjuntos, pueden
>>> contener información privilegiada y/o confidencial que está dirigida
>>> exclusivamente a su destinatario. Si usted recibe este mensaje y no es el
>>> destinatario indicado, o el empleado encargado de su entrega a dicha
>>> persona, por favor, notifíquelo inmediatamente y remita el mensaje original
>>> a la dirección de correo electrónico indicada. Cualquier copia, uso o
>>> distribución no autorizados de esta comunicación queda estrictamente
>>> prohibida.
>>>
>>> CONFIDENTIALITY NOTICE: This e-mail and the accompanying document(s) may
>>> contain confidential information which is privileged and intended only for
>>> the individual or entity to whom they are addressed. If you are not the
>>> intended recipient, you are hereby notified that any disclosure, copying,
>>> distribution or use of this e-mail and/or accompanying document(s) is
>>> strictly prohibited. If you have received this e-mail in error, please
>>> immediately notify the sender at the above e-mail address.
>>>
>>
>>
>>
>