Re: [eigen] Re: efficient trace of product of symmetric matrices

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

On Mon, Feb 20, 2012 at 11:28 PM, Patrick <pat05.mlst@xxxxxxxxx> wrote:
> Hi everyone,
>
> I'm wondering whether it is currently possible to compute
>
> tr(A*B),   where   A \in Sym(n,n), B \in Sym(n,n)
>
> efficiently using the development branch of Eigen3.
> Since A and B are symmetric and
>
> tr(A*B) = sum_ij a_ij * b_ij,
>
> the best way I can think of doing this is
>
> double tr=0;
> for(size_t i=0; i<A.rows(); ++i)
> {
>     tr += A(i,i) * B(i,i);
>     for(size_t j=0; j<i; ++j)
>         tr += 2 * A(i,j)*B(i,j);
> }
>
> I noticed that the triangularView methods do not expose the array
> interface, so am I correct in assuming I cannot vectorise this using
> triangular views?

what you need here is a lazy triangular * triangular product, and
indeed there is no such thing yet. Nevertheless you can still write
your inner loop using an Eigen expression and get vectorization if A
is row major and B column major.

> Furthermore, when trying to make a matrix symmetric by calling
>
> A.triangularView<StrictlyUpper>() = A.triangularView<StrictlyLower>(),
>
> I noticed that no assignment was taking place, i.e. the matrix A remained as
> is.
> Is this the expected behaviour?

Yes that's the intended behavior. You have to transpose the right hand side:

A.triangularView<StrictlyUpper>() =
A.triangularView<StrictlyLower>().transpose(),

Actually, you can even do:

A.triangularView<StrictlyUpper>() = A.transpose();

The triangularView on the lhs behave like a writing mask.

cheers,
gael.

>
> Thanks.
>
> Kind regards,
>
> Patrick