[eigen] Re: efficient trace of product of symmetric matrices
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: [eigen] Re: efficient trace of product of symmetric matrices
• From: Patrick <pat05.mlst@xxxxxxxxx>
• Date: Mon, 20 Feb 2012 22:28:53 +0000
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; bh=4UNYvwlhG4h75AN9IiW19PGcYdTrYCicVRy9mDm87Eg=; b=j+M4pYAzJxYF6Dgu3p6NQto88ekurEwqGLWHMOG23mVoBuPkxXOL+cfTO/mZNklvHI YCP1k0HzZ7QmSaV7Pl7Fq302Sv/v7EEvH55n6RffkrGPXfmUSAlAr4syGiKoH7yffP+1 3vGk32IrrfnxOWj6kBmTC4T/nvZS26zkQ1O3c=

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?

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?

Thanks.

Kind regards,

Patrick

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