Re: [eigen] part<SelfAdjoint> in Eigen3

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

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