Re: [eigen] Allocation policy of eigen decompositions

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

2010/4/21 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Tue, 20 Apr 2010, Benoit Jacob wrote:
>
>> 2010/4/20 Adolfo Rodríguez Tsouroukdissian
>> <adolfo.rodriguez@xxxxxxxxxxxxxxxx>:
>>>
>>> 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.
>
> I'm definitely interested. However, I didn't have time to reply earlier. I
> don't see a problem with the changes proposed by Benoit. I'm glad to hear
> that the template Options parameter in JacobiSVD is on the way out.
> ReturnByValue is your friend and easy to use.
>
>> 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().
>
> I guess I don't understand how HouseholderSequence works. In fact, I thought
> about using it when rewriting HessenbergDecomposition, but I refrained from
> doing so because I had trouble understanding the code.

A HouseholderSequence is just a wrapper around an existing matrix
where the essential parts of the householder vectors are stored in
some triangle, and around a separate vector storing the "householder
coefficients". Thus creating a householderSequence object doesn't
cause any malloc, just like creating a Block expression or a
TriangularView doesn't cause any malloc.

> In this code,
>
>   MatrixXd A = MatrixXd::Random(100,100);
>   HouseholderQR<MatrixXd> houseOfA(A);
>   MatrixXd Q = houseOfA.householderQ();
>
> I'd have thought that the last line constructs a temporary
> HouseholderSequence objects

Yes it does

> which holds a m_vectors member (matrix of size
> 100 x 100),

It only holds a _reference_ to such a matrix.

Here's how the source code looks like (HouseholderSequence.h) :

template<typename VectorsType, typename CoeffsType, int Side> class
HouseholderSequence
{
  ...

  protected:
    typename VectorsType::Nested m_vectors;
    typename CoeffsType::Nested m_coeffs;
    bool m_trans;
    int m_actualVectors;
    int m_shift;
};

See the ::Nested here. typename VectorsType::Nested is equivalent to
ei_nested<VectorsType>::type which in the case of a plain matrix type,
will give you a reference type. See XprHelper.h.

> calls its evalTo() method, and throws the temporary object away.
> Isn't there a malloc needed to allocate the memory in which the m_vectors
> member is stored?

No malloc, since m_vectors is just a reference.


>
>
> And here is my list of issues that may require API changes.
>
>
> Chaining decompositions
> =======================
>
> As Adolfo says, ComplexSchur uses HessenbergDecomposition. If you call
> ComplexSchur::compute(matrix), then the ComplexSchur object creates a
> HessenbergDecomposition object, calls its compute() method, and copies the
> matrices computed into its own member variables. One problem is that
> creating the HessenbergDecomposition object causes a memory allocation. That
> is the issue that Adolfo is concerned with.

Once HessenbergDecomposition's matrixQ() method just returns a
householder sequence, and once  HessenbergDecomposition gets the
pre-allocation API from Adolfo's patch, this is a problem solved,
right?

> Another problem is that this uses more memory than necessary. The
> HessenbergDecomposition object needs memory and the ComplexSchur object
> needs memory. Ideally, they should use the same memory, instead of copying
> the matrix. What is the best way to do that? Perhaps ComplexSchur can use
> the memory under control of HessenbergDecomposition by casting away the
> const-ness of the const reference returned by
> HessenbergDecomposition::packedMatrix(), but this is not very kosher.

We could just add a non-const-qualified getter method, packedMatrix().
No big deal. If people want to ensure that their HessenbergDecomp is
treated as a constant object, they just have to mark it const, which
will ensure that such a non-const getter doesn't get called.

> The same issue arises in other contexts in the Eigenvalues module: RealSchur
> uses HessenbergDecomposition, EigenSolver uses RealSchur, ComplexEigenSolver
> uses ComplexSchur. In all these cases there is potential for re-using
> memory, reducing the memory footprint. Whether it's worthwhile to save
> O(n^2) memory in an O(n^3) algorithm is another matter ... perhaps not the
> most pressing issue.

Not pressing indeed since, if you agree with my solution, this is only
a matter of adding a non-constant getter, which doesn't require an API
break.


> Merging EigenSolver and ComplexEigenSolver
> ==========================================
>
> My other big question is whether to merge the EigenSolver and
> ComplexEigenSolver classes.

I am totally incompetent here so will leave it to others...

That doesn't sound trivial as they're doing different things,
pseudoeigenvectors versus actual-eigenvectors, etc. As you mention.

Benoit

> Intuitively, I should say yes, and that's also
> what the ToDo list on the wiki says. However, the interface of EigenSolver
> and ComplexEigenSolver are rather different. EigenSolver stores the
> eigenvalues and the pseudo-eigenvectors, so EigenSolver::eigenvectors()
> returns MatrixType (this can be replaced via the ReturnByValue mechanism but
> that does not matter for this discussion). ComplexEigenSolver stores the
> eigenvalues and the eigenvectors, so ComplexEigenSolver::eigenvectors()
> returns a const ref to MatrixType. How do we merge them?
>
> One possibility is to have eigenvectors() in the merged EigenSolver class
> return a MatrixType if the scalar type is real and a const ref to MatrixType
> if the scalar type is complex. This is similar to eval() returning sometimes
> a matrix and sometimes a const ref. However, it may be confusing to users..
>
> The second possibility is to have the merged EigenSolver class store the
> eigenvectors and not the pseudo-eigenvectors in the real case. However, this
> costs twice as much memory because the eigenvectors are complex and the
> pseudo-eigenvectors are real. As the real case is the important case, this
> is not ideal.
>
> Then there is the minor issue that (real) EigenSolver has a
> pseudoEigenvector() method and ComplexEigenSolver has not. However, we can
> include a pseudoEigenvector() method in the merged EigenSolver class, and
> have it return the eigenvectors in the complex case. I don't think this is a
> problem.
>
> There is a third possibility: not to merge (real) EigenSolver and
> ComplexEigenSolver.
>
> For bonus points, discuss merging RealSchur and ComplexSchur.
>
>
> Cheers,
> Jitse
>
>
> PS: Adolfo, the Eigenvalues part of your patch looks good, at first sight..
> Incidentally, how do you test that there are no malloc's?