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*: "ESCANDE Adrien 222264" <adrien.escande@xxxxxx>*Date*: Thu, 29 Jul 2010 17:14:08 +0200*Thread-index*: AcsvHxib6WrRC3CxR72mbTSPawvw0AADTsDg*Thread-topic*: [eigen] part<SelfAdjoint> in Eigen3

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**:**[eigen] part<SelfAdjoint> in Eigen3***From:*ESCANDE Adrien 222264

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Do we need geometry refactoring?** - Next by Date:
**[eigen] Feature request: TriangularView::conjugate()** - Previous by thread:
**Re: [eigen] part<SelfAdjoint> in Eigen3** - Next by thread:
**[eigen] Feature request: TriangularView::conjugate()**

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