Re: [eigen] Allocation policy of eigen decompositions |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Allocation policy of eigen decompositions*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Wed, 21 Apr 2010 12:46:07 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:received:message-id:subject:from:to:content-type :content-transfer-encoding; bh=PXYNAR74Mtq8TkwmQiPQ8sO3x/+apgoDiC+ykos6/Yg=; b=jHQmdzjbg5BxTPJ1beLLJdSwxMbCt9Z+mmDqOJZirqwYaGEoGOmJ1wokDB1qPuRzb0 /8i/ZC8NVGnMEmMR1FZfzHjfyeidu5c2WzmfdlLkRk1rnOT/B3RZUXHCu76t3keQjgoN Hpfbc1JV53aEuMGBg553040wdX9+sQJxGA+EE=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=utCcFhb7w10tsxOxCdJLFX0Q1zevfvDYVShq/SGzKKmJXSNTNJENoklCTa7vpB0KeK hZcilCQZM+5vYrCrUQPSgCeuTqbqaWwqky1dcwlWegc28Wu+eZwPvioRPRVzrLBUWW/s QCFr6GuWuTNk19AK387HlIfu20mkVToYS1yyY=

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?

**References**:**Re: [eigen] Allocation policy of eigen decompositions***From:*Benoit Jacob

**Re: [eigen] Allocation policy of eigen decompositions***From:*Adolfo Rodríguez Tsouroukdissian

**Re: [eigen] Allocation policy of eigen decompositions***From:*Benoit Jacob

**Re: [eigen] Allocation policy of eigen decompositions***From:*Adolfo Rodríguez Tsouroukdissian

**Re: [eigen] Allocation policy of eigen decompositions***From:*Benoit Jacob

**Re: [eigen] Allocation policy of eigen decompositions***From:*Adolfo Rodríguez Tsouroukdissian

**Re: [eigen] Allocation policy of eigen decompositions***From:*Adolfo Rodríguez Tsouroukdissian

**Re: [eigen] Allocation policy of eigen decompositions***From:*Benoit Jacob

**Re: [eigen] Allocation policy of eigen decompositions***From:*Jitse Niesen

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Compilation errors introduced by recent changes. (aligned stl types & MatrixType rename)** - Next by Date:
**Re: [eigen] Allocation policy of eigen decompositions** - Previous by thread:
**Re: [eigen] Allocation policy of eigen decompositions** - Next by thread:
**Re: [eigen] Allocation policy of eigen decompositions**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |