Re: [eigen] part<SelfAdjoint> in Eigen3 |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] part<SelfAdjoint> in Eigen3*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Mon, 16 Aug 2010 09:00:10 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:in-reply-to :references:date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=sQauhzh2y7w/1D77revmVZLt58rWff7IC1WCymgV6uY=; b=fZ3SauGMUv/5HJNKdklmDu2ek/Q8epiZ04zFunBzFv69LqZ3N+cUOOBFQrOrZqOvOE 4DyTj0AdavRoNXGdhkhRnGFMGUsLVEWjs3dwRlW2L3x0jeRJLb1weC5ixsFpjaZNJD6o bt7jyXcQT5Uqjz3tixmskMIUx+rQ2QZpRNwJE=*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=kyKqS7iV8OqpxwAU3LngwnObaUEE1KpM4LSxzf3kyEQHA8eTbY4wdZDg1XKZ6I5Xz4 EtVz+428kwmiiUBorWW4boqxbWn/0xdwG6+t/dIZWqcCQZuiYi22NFiTD/MbxFLKoIxy 1eJeF92CIw53mDP8umACGk1dc8r4mjaLSoJIg=

2010/8/16 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>: > The selfadjointView mechanism aims at viewing a triangular part > (either upper or lower) as a selfadjoint matrix. Overall Eigen > (products, decompositions, solvers) you can deal with such selfadjoint > matrices where only one half is meaningful. So there is very little > use to add methods to fill both triangular parts at once. If you > really want to have the full matrix you can still do: > > // fill only one half, e.g.: > m.triangularView<Upper>() = v + v.adjoint(); > m.selfadjointView<Upper>().rankUpdate(u); > > > // get the full matrix: > m.triangularView<Lower>() = m.adjoint(); > > > Finally, for large enough matrices, we could easily detect at runtime > products like: > > m = u.adjoint() * u > > and automatically optimize it, just like ATLAS's xGEMM routines do. OK to do that for dynamic-size matrices, but that still leaves out the case of fixed-size matrices. However, since I don't think we have any optimize "rankUpdate" code for fixed-size matrices, I guess that doesn't matter for now. We can do the optimization you describe now for dyn-size and worry later about fixed-size (consider adding a method like I suggested). Benoit > > cheers, > > gael > > > > On Thu, Jul 29, 2010 at 5:14 PM, ESCANDE Adrien 222264 > <adrien.escande@xxxxxx> wrote: >> The problem is that there doesn't seems to be an Eigen3 way to indicate that the result of an operation is a selfadjoint matrix and thus that the computation can be optimized (i.e. only about half of the coefficients need to be really evaluated). >> n.selfadjointView<UpLo>() does not have an operator= (or so my compiler tells me when I try to write Eq (1)). >> No way thus to optimize (1) or something like n = m.transpose()*weight..asDiagonal()*m, with weight a vector. (No direct way at least, I could explicitly compute a triangular part of n then copy it into the other). >> >> Template parameter UpLo tells to consider only the coefficients in the upper/lower part of the matrix and deduce the other coefficients from them as can be seen in this example. >> >> MatrixXd M = MatrixXd::Random(8,8); >> M.triangularView<StrictlyUpper>().setZero(); >> std::cout << M << std::endl << std::endl; >> //std::cout << M.selfadjointView<Upper>() << std::endl; //FAILS >> //MatrixXd M2 = M.selfadjointView<Upper>(); //FAILS >> MatrixXd M2 = M.selfadjointView<Lower>()*MatrixXd::Identity(8,8); >> std::cout << M2 << std::endl << std::endl; >> >> The example underlines two problems: you can print out a View (triangular or selfadjoint), and while you can assign a TriangularView to a Matrix, you can't assign a SelfAdjointView. >> >> >> Adrien >> >> >> >> -----Message d'origine----- >> De : Listengine [mailto:listengine@xxxxxxxxxxxxxxxxx] De la part de Benoit Jacob >> Envoyé : jeudi 29 juillet 2010 15:08 >> À : eigen@xxxxxxxxxxxxxxxxxxx >> Objet : Re: [eigen] part<SelfAdjoint> in Eigen3 >> >> 2010/7/29 ESCANDE Adrien 222264 <adrien.escande@xxxxxx>: >>> Hello, >>> >>> >>> >>> Eigen2 documentation mentions the way to obtain an optimized computation for >>> an expression evaluating to a selfadjoint matrix: >>> >>> n.part<SelfAdjoint>() = m+m.adjoint ; (1) >>> >>> n.part<SelfAdjoint>() = (m*m.adjoint()).lazy(); (2) >>> >>> >>> >>> I didn't find the way to have that in Eigen3 (MatrixBase>Derived>::part is >>> still there but flagged as deprecated). The Porting from Eigen2 to Eigen3 >>> page does not mention this case (it gives only the translation for >>> part<SelfAdjoint|Upper> and part<SelfAdjoint|Lower>). Is there a new direct >>> way to perform such computations? >> >> The new method is selfadjointView(). >> >> n.selfadjointView<Lower>() >> n.selfadjointView<Upper>() >> >> >>> >>> And what would be the best way to write (2) when m itself is selfadjoint? >>> >> >> Gael would know better, but I think that your best bet is rankUpdate() here. >> >> n.setZero(); >> n.selfadjointView<Upper>().rankUpdate(m); >> >> I can't see any way to do that using arithmetic operators, but that >> seems OK as indeed the API for doing that should only take one operand >> m. >> >> What I am a bit more puzzled about is that this API forces the user to >> pass Upper... I need an explanation as to why this is useful ? >> >> Maybe it would be useful to add a productByAdjoint(), perhaps >> abbreviated as .xxt(), method in MatrixBase so the user could do: >> >> n = m.productByAdjoint(); >> >> ? >> >> Benoit >> >>> >>> >>> Adrien >>> >>> >>> >>> >> >> >> >> >> > > >

**References**:**Re: [eigen] part<SelfAdjoint> in Eigen3***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] InnerPanel change mis-detects alignment?** - Next by Date:
**Re: [eigen] [eigen3] .cwise() vs .array() and static assertion** - Previous by thread:
**Re: [eigen] part<SelfAdjoint> in Eigen3** - Next by thread:
**Re: [eigen] Expression templates**

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