Hi guys,

is there a way to extract the symmetric part of a squared matrix more efficiently than M1=(M+M.transpose())/2.0 ?

Thanks for you help!

sorry but I don't really understand what you wanna do. Could you be more precise?

Hi, sure, I'll try to be more precise.

I'm running some code that outputs a tensor T which I represent as
Eigen::Matrix<double,dim,dim> T

In continuum mechanics you often need to decompose a tensor in its symmetric and antisymmetric parts. In this case I'm interested in the symmetric part of T, which I calculate simply as
Ts=(T+T.transpose())/2.0

Since I need to do this operation many times I was wondering if there is a more efficient way to extract the symmetric part of my tensor T.

Thank you.

Note that there is a SelfAdjointView. I know this is not exactly the same thing, but you have real matrices, so self-adjoint is the same as symmetric.

Could an expert confirm if this would work?

I think that here n2k is referring to the fact any matrix (or tensor) can be represented as a sum of symmetric and anti-symmetric matrices:

A = (A+A^T)/2 + (A-A^T)/2

A compact way to achieve and store this is to use the triangularView, e.g.:

B.diagonal() = A.diagonal();
B.triangularView<StrictlyLower>() = (A.triangularView<StrictlyLower>() + A.triangularView<StrictlyUpper>().transpose())*0.5;
B.triangularView<StrictlyUpper>() = (A.triangularView<StrictlyUpper>() - A.triangularView<StrictlyLower>().transpose())*0.5;

Then you can get an expression of (A+A^T)/2 as B.selfadjointView<Lower>(), but unfortunately there is no such anti-symmetric view yet. So you would have to do things by hand using B.triangularView<Upper>() and/or -B.triangularView<Lower>().

If you tell us more about you wanna do with such decomposition we could give you more relevant hints maybe...

ggael wrote:I think that here n2k is referring to the fact any matrix (or tensor) can be represented as a sum of symmetric and anti-symmetric matrices:

A = (A+A^T)/2 + (A-A^T)/2

This is exactly what I'm talking about.

ggael wrote:If you tell us more about you wanna do with such decomposition we could give you more relevant hints maybe...

ok I'll try. I'm writing a finite element library for defect mechanics in crystals that uses Eigen for linear algebra (thanks for the great job btw). Speed is crucial for me and I spend most of the time integrating (quadrature) some tensor, let's say A, which I represent as:

Code: Select all

Eigen::Matrix<double,dim,dim> A;


It turns out that A is a symmetric tensor and can be written as

Code: Select all

A = 1/2*integral (  LongExpression+ LongExpression^T )


Taking advantage of linearity I figured that the following would be faster

Code: Select all

A1 = integral (  LongExpression )
A = 1/2*(A1+A1^T)=symmetricPart(A1)


so I end up extracting the symmetric part many times, which is why I asked if there is an efficient way to do it. I guess an alternative (which I'm going to try) would be

Code: Select all

A.part<Eigen::SelfAdjoint>()= integral (  LongExpression+ LongExpression^T )


... but I'm not that familiar with part/SelfAdjoint.

Thanks again!

Ah, ok so your are simply looking for an efficient way to compute (A + A^T)*0.5. So what I wrote in my previous message was for the devel branch only. With the 2.0 branch you can do:

A.part<SelAdjoint>() = (A1 + A1.transpose())*0.5;

This will work and set the complete matrix A while computing only one triangular part.

But now what do you want to do with A ? Maybe you don't need to store both triangular part of the matrix A. Actually here the the random write access to write to symmetric are likely to be quite costly... For instance you could even compute A in place like this:

A1.part<StrictlyLower>() = (A1.part<StrictlyLower>() + A1.transpose().part<StrictlyLower>())*0.5;

which should be more efficient.