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

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

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