Re: [eigen] Allocation policy of eigen decompositions

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

2010/4/20 Adolfo Rodríguez Tsouroukdissian <adolfo.rodriguez@xxxxxxxxxxxxxxxxx>:
> I wanted to bring up some issues concerning the Eigenvalues module.

I am adding Jitse in CC to ping him, although he is reading this list
anyway; he has been improving the Eigenvalues module a lot recently so
might be interested.

> All 7
> decompositions now have a problem size constructor, but the final goal,
> which is to prevent heap allocations in the compute() method when using
> (smallish) dynamic matrices + the problem size constructor, seems a bit more
> difficult to achieve than with the other modules. These are the main
> holdbacks:
>
> 1. Methods returning matrices by value
>
> Most other eigen decompositions have output getters that return
> reference-to-const matrices/vectors, which is very convenient. In the cases
> where an output is currently not a class member, but constructed on demand
> (as in HessenbergDecomposition<MatrixType>::matrixQ() ), a template Options
> parameter à la JacobiSVD would do the trick and enable returning a
> ref-to-const, n'est-ce pas?.

But we are moving away from this approach, because template parameters
mean larger binary code when using more than one variant, and make for
rather complex APIs. I'm going to do that change in JacobiSVD Real
Soon Now (tm).

Plainly returning by value is also something we'll move away from. The
direct improvement on that is returning a ReturnByValue proxy object,
i.e. a light object describing the operation to be done, letting the
actual computation occur when it's assigned to a destination matrix.
Since such a destination matrix is created by the user, the user is in
control of the allocation, which should be good from your perspective,
right ? :)

For an example of ReturnByValue, see for example
PartialPivLU::solve(). Notice that ei_solve_retval inherits
ReturnByValue (see in Eigen/src/misc/Solve.h).

That said, in the case of
HessenbergDecomposition<MatrixType>::matrixQ(), we can and should do
much better than that: return a HouseholderSequence object, describing
the set of householder vectors stored in the matrix. That is not only
much more efficient algorithmically, that also avoids the malloc. For
example, see HouseholderQR::householderQ().

Finally, if eventually some decomposition class really had to allocate
a matrix for the result or an intermediate result of an operation, one
could make that matrix a member of the class... it's just that this
should (obviously) be a last-resort solution.

> Consider the particular example of
> ComplexSchur<MatrixType>::compute(), which uses internally a Hessenberg
> decomposition. Now I've made HessenbergDecomposition<MatrixType>::compute()
> heap-allocation-free, but then as soon as I invoke
> HessenbergDecomposition<MatrixType>::matrixQ() I have an unwanted (and
> preventable) heap allocation.

Yep, that'll go away as soon as this matrixQ() method is replaced by a
householderQ() method as explained above :)

> As a sidenote, there is a comment in line 170 of Tridiagonalization.h that
> states:
>
> // FIXME should this function (and other similar ones) rather take a matrix
> as argument and fill it ? (to avoid temporaries)
>
> but from the perspective of maintaining consistency with the rest of Eigen,
> maybe returning ref-to-consts is best (?)

This must be an old comment, from the days before we had proper Band
matrices in Eigen. Tridiagonalization needs to be rewritten to imitate
UpperBidiagonalization, which is modern and uses the BandMatrix
framework. This matrixT() method then just returns a reference to the
stored BandMatrix.

> 2. Prevent the creation of temporaries
>
> Although the simpler cases have already been addressed, there are still a
> few less-trivial cases where it's better to ask for an opinion. For example,
> in SelfAdjointEigenSolver<MatrixType>::compute(...) there is this line:
>
> TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag,
> computeEigenvectors);
>
> Which internally allocates a Tridiagonalization solver. What would be --in
> your opinion-- the best workaround in terms of API consistency and code
> duplication prevention?.
>
>
> That's it for now. Of the two issues described above, 1. seems to be the
> more pressing one, since it potentially involves API changes. As always,
> thanks in advance for your advice.

Ok so I'm really tired after having spend 5 hours fixing the bug that
you reportes (thanks by the way) so if you tell me that 2. can wait a
bit, I'll wait a bit :)

Benoit

>
> Adolfo
>
>
>
>
> 2010/4/20 Adolfo Rodríguez Tsouroukdissian <dofo79@xxxxxxxxx>
>>
>>
>> On Mon, Apr 19, 2010 at 7:36 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
>> wrote:
>>>
>>> Sorry, my bad! I added this eval() for debugging and forgot to remove it.
>>>
>>> Can you retry now? This should remove 1 heap allocation. I really
>>> don't know where the 3 others can be coming from, but that would be a
>>> bug that we (I...) introduced, so you don't have to block your stuff
>>> because of it.
>>
>> Removing eval() did the trick. I was seeing four allocations because the
>> call to applyHouseholderOnTheLeft is inside a loop :s, so problem solved..
>>
>> Adolfo
>>
>>>
>>> Benoit
>>>
>>>
>>> 2010/4/19 Adolfo Rodríguez Tsouroukdissian
>>> <adolfo.rodriguez@xxxxxxxxxxxxxxxx>:
>>> > Some updates,
>>> >
>>> > The patch is coming out quite well, but there is something that I'd
>>> > like to
>>> > ask the devs about firs. Up until this morning I had all decompositions
>>> > except those in the Eigenvalues module (will talk about that in a
>>> > separate
>>> > email) working as expected, then I pulled in recent changes
>>> > (storageorders
>>> > branch included) and something changed in the behavior of the
>>> > applyHouseholderOnTheLeft method of MatrixBase<Derived>.
>>> >
>>> > Line 101 of Householder/HouseHolder.h states:
>>> >
>>> > tmp.noalias() = essential.adjoint().eval() * bottom;
>>> >
>>> > ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
>>> > When using dynamic matrices, the highlighted part performs four heap
>>> > allocations, and since all QR decompositions use
>>> > applyHouseholderOnTheLeft,
>>> > their compute() method ends up heap-allocating as well, which is what
>>> > I'm
>>> > trying to avoid. My question is then, can these allocations be
>>> > prevented?,
>>> > and if so what would be the best way?. I would also appreciate some
>>> > insight
>>> > into why these heap allocations weren't occurring before.
>>> >
>>> > Thanks in advance,
>>> >
>>> > Adolfo
>>> >
>>> >
>>> > --
>>> > 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.
>>> >
>>>
>>>
>