Re: [eigen] SelfAdjointView<> as LHS expression for optimizations?

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

For the implementation, see what I wrote above, it is just a matter of
a few cosmetic changes in internal::selfadjoint_product.... so don't
worry !

gael

On Tue, Nov 9, 2010 at 1:18 AM, Jose Luis Blanco
<joseluisblancoc@xxxxxxxxx> wrote:
> Hi Gael, Benoit,
>
> Thanks for the pointers!
>
> Giving it a second thought, I also think a RHS approach would be
> clearer and perhaps also easier to implement.
>
> I like the:
>  result  =  (H * C * H.transpose()).assumeSymmetric();
>
> method, since it would be more generic than the HCH^t operation.
>
> But on twist() as a method, it would probably be the easiest way to
> go, although the operators I know as "twists" are for converting from
> SO(3) to so(3) and such, and I don't see the resemblance at first
> sight...
>
> If you give me some basic starting points to start with, I could *try*
> implementing such "C.twist(H)" method.
> Naively I would go and write it as part of MatrixBase<> as something
> like http://reference.mrpt.org/0.9.1/ops__matrices_8h-source.html#l00469
> , but obviously I can do it much better with Eigen, so please let me
> know another "similar" methods already in Eigen and/or docs I can
> study to implement it more efficiently taking advance of
> vectorization...
>
> In the topic I work at, a very common operation in inner loops is
> H*C*H^t with C a cov. matrix and H a Jacobian, that's the why of my
> interest in optimizing the operation as much as possible.
>
> JL
>
>
> On Tue, Nov 9, 2010 at 12:07 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> 2010/11/8 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>> On Fri, Nov 5, 2010 at 8:49 PM, Jose Luis Blanco
>>> <joseluisblancoc@xxxxxxxxx> wrote:
>>>> Hi all,
>>>>
>>>> The operation I have in mind is:
>>>>  RESULT = H * C * H.transpose()
>>>> with C symmetric and squared.
>>>>
>>>> The point is that two "problems" (well... I have to admit, small
>>>> ones!) occur with a naive Eigen expresion "H * C * H.transpose()":
>>>>
>>>> 1) Many multiplications could be saved by realizing that the result is
>>>> exactly symmetric. I know the implications may be far reaching, but,
>>>> could SelfAdjointView<> be extended so it can be used to store
>>>> results? I mean, something like:
>>>>
>>>> Eigen::MatrixXd RESULT;
>>>> RESULT.selfadjointView<Eigen::Lower>() =
>>>> H*C.selfadjointView<Eigen::Lower>()*H.transpose();
>>>>
>>>> or an alternative that wouldn't require LHS  = operator for SelfAdjointView<>:
>>>>
>>>> Eigen::MatrixXd RESULT =
>>>> (H*C.selfadjointView<Eigen::Lower>()*H.transpose()).selfAdjointView<Eigen::Lower>();
>>>>
>>>> Any ideas about the difficulty of implementing any of these two ideas?
>>>> Or... is there anything similar already done?
>>
>> I would say that a really good API here would only be about the
>> right-hand side, since the left-hand side is irrelevant to the fact
>> that one can do this optimization. We could either (that's the
>> simplest to implement) have a specialized method for that, which
>> according to mathematicians jargon could be called twist(), like
>>
>>   result = C.twist(H)
>>
>> Or, once bug 99 is implemented, we could modify the product expression
>> after it's constructed to tell Eigen that we know it's symmetric:
>>
>>   result  =  (H * C * H.transpose()).assumeSymmetric();
>>
>> In any case I would prefer to try solutions involving only the RHS
>> before we look at a LHS API.
>>
>> Benoit
>>
>>>
>>> Good question. Actually I don't any particular optimizations beside
>>> computing only half of the second product. Basically, this would be
>>> implemented as follow:
>>>
>>> Temp = H * C.selfadjointView<Lower>();
>>> Result.triangularView<Lower> = Temp * H.adjoint();
>>>
>>> The second line currently does not work, however this is very easy to
>>> do by extending our internal selfadjoint product routine to compute
>>> half of a general A*B product instead of handling AA* or A*A only.
>>> Actually this routines cannot exploit the fact that both sides come
>>> from the same matrix anyway....
>>>
>>> Regarding the API we could simply let the user write:
>>>
>>> Result.triangularView<Lower> = H * C.selfadjointView<Lower>() * H.adjoint();
>>>
>>> For consistency reason,
>>>
>>> Result.selfadjointView<Lower> = H * C.selfadjointView<Lower>() * H.adjoint();
>>>
>>> would do the same.
>>>
>>> gael
>>>
>>>>
>>>> 2) Apart from efficiency, there is this numerical accuracy problem of
>>>> the resulting matrix not being *exactly* symmetric... a*b*c != c*b*a
>>>> for most architectures. Ok, the errors are tiny, but it may lead to
>>>> inconsistencies that some time ago affected the Cholesky
>>>> implementation I used to use... that's why I worry about all this in
>>>> the first place!
>>>>
>>>>
>>>> Best,
>>>> JL
>>>>
>
>
>