RE: [eigen] part<SelfAdjoint> in Eigen3

["ESCANDE Adrien 222264" <adrien.escande@xxxxxx>]

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