Re: [eigen] Re: efficient trace of product of symmetric matrices
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] Re: efficient trace of product of symmetric matrices
• From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
• Date: Tue, 21 Feb 2012 23:20:16 +0100
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; bh=wXxI7nzY+UC6mAKH7Vmu96fsppUkoEZHo1pNnG16hv8=; b=D19+onOun4VIJHY5VU8uXoDP5HOQl1CBf2Y7t5FIKJ87Hmjt99eIUfbGJB0/ebDvnI HDOSZdUJw8sEHMrNVGp82+rRc5I98pHV+uVMG9nethVXWoMlVAkvuuXEqC9zQ5rL5fb4 hnfoH5XLKeJZmrfxyRZ6Nix4uozUQ2QTylBhw=

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



• References:

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/