| Re: [eigen] Allocation policy of eigen decompositions |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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.
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 which holds a m_vectors member (matrix of size
100 x 100), 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?
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.
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.
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.
Merging EigenSolver and ComplexEigenSolver
==========================================
My other big question is whether to merge the EigenSolver and
ComplexEigenSolver classes. 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?