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*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Mon, 16 Aug 2010 10:12:09 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type :content-transfer-encoding; bh=9G7DKD0lGlIczPV/a1nPFLcEfpovMufAukaM6mqVWV0=; b=SLgLb1F1cfnsdeS7acSN4d8XcbPF03mDYYpmxm5wHRu5LW+M8jypLz6xU99trfFTmu 5PqqkLg39iRgilhQ+w/lsZ0748MfTphvXp/hleqC4rrOYgNfCVe4qbblqj8rTPwPGeQn amNy0ejiCynOXtSw7l+v9HTg8FCsos1aShwCU=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; b=C+vvUx2IBWHJWindSRMkPLTSWpmBbJjOsxZjIpkP3Dl0pV2VxG5jLm9aRiW8Eh4vIm Nib+s7yyxNmZQ89AKrkyGyjKA9Lnp7i3fEEWJ/6aXI/y5Gl2APZBs1Z6JNbdysxVulJg z8MypSfNLDmKfQwDgtbFrwtteCQxrULXCa5lo=

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. 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 >> >> >> >> > > > > >

**Follow-Ups**:**Re: [eigen] part<SelfAdjoint> in Eigen3***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] InnerPanel change mis-detects alignment?** - Next by Date:
**Re: [eigen] Expression templates** - Previous by thread:
**Re: [eigen] Matrix product crashes when compiled with MSVC 2010 in release** - Next by thread:
**Re: [eigen] part<SelfAdjoint> in Eigen3**

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