Re: [eigen] part<SelfAdjoint> in Eigen3

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

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