|Re: [eigen] SelfAdjointView<> as LHS expression for optimizations?|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] SelfAdjointView<> as LHS expression for optimizations?
- From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Tue, 9 Nov 2010 09:20:09 +0100
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type :content-transfer-encoding; bh=Zygdk00TwgNKGRstZ6GAYIoOS/BK64aWGeH1RgF4NBc=; b=jOggIPqnUM4fhGIg8/+HNhdNctIVGW7QNRhLlrjLnsavFZ0wMy9xpvWx84WB79TAJs 20dK9fVR5a4GraG9h4EuQEAP4uG9wSvM1PV/bT+M+UDMSurSX1pRBSkeJJh83qfP2mge 4lRqn8eaE0vxwhXMbBMGA5JfZdw0dUbzrjeWY=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; b=xTRQUgsFh8C/jkKRb276SbaIGXDhtPhW/r8EVQe7t8RqTeQVEU+fti3s1/jzNZYjTJ 7tJCCS7n8Q0uWFZXWOL1i/n+l+gLd1QmmnJTGkn9KkSLgCnpyg+Fbm3M9neQfZv0OxHM WpfzCXMsPCo/ZBzcbO9rdmzjkEoI7GW4KVXYg=
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
On Tue, Nov 9, 2010 at 1:18 AM, Jose Luis Blanco
> 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
> 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
> 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.
> 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>() =
>>>> or an alternative that wouldn't require LHS = operator for SelfAdjointView<>:
>>>> Eigen::MatrixXd RESULT =
>>>> 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.
>>> 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.
>>>> 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!