Exploiting symmetry of A * S * trans(A)

[Mikhail Katliar]

The result of A * S * A^T is a symmetric matrix provided that S is symmetric. However, in the following code Blaze evaluates the right-hand side of the minus-assignment into a non-symmetric temporary result:

using Matrix = blaze::DynamicMatrix<Real>;
using SymmetricMatrix = blaze::SymmetricMatrix<Matrix>;

SymmetricMatrix stateCovariance_;
SymmetricMatrix S_;
Matrix K_;

// ...

stateCovariance_ -= K_ * S_ * trans(K_); // rhs is evaluated to a non-symmetric temporary matrix, although it is always symmetric

The expression of form A * S * A^T are quite common in linear algebra. Is it possible to handle them in Blaze in such a way that the symmetry is exploited?

[Klaus Iglberger]

Hi Mikhail!

I have good news for you: The feature already exists in the form of the declsym() function:

stateCovariance_ -= declsym( K_ * S_ * trans(K_) );

For more information, see the wiki or issue #35.

Best regards,

Klaus!

[Klaus Iglberger]

Duplicate of #35.

[Mikhail Katliar]

Hello Klaus,

Thank you for the good news! But I thought that declsym() is meant to be used when it is not possible to determine that the expression is symmetric at compile-time. In case of, for example, addition of two symmetric matrices, Blaze automatically determines that the result is symmetric and declsym() is not necessary. Similarly to this case, if S is a SymmetricMatrix then A*S*A^T is also always symmetric. This is known at compile time and it looks to me that it could work without declsym(). Or do I misunderstand something?

Best regards, Mikhail

[Klaus Iglberger]

Hi Mikhail!

Unfortunately it is not in general possible to determine that a given matrix multiplication will result in a symmetric matrix at compile time. It is not the type of the matrices but the actual instances that provide information about symmetry. Consider the following example:

using Matrix = blaze::DynamicMatrix<Real>;
using SymmetricMatrix = blaze::SymmetricMatrix<Matrix>;

SymmetricMatrix S_;
Matrix K_;
Matrix L_;

using ExprType1 = decltype( K_ * trans(K_) );
using ExprType2 = decltype( K_ * trans(L_) );
static_assert( std::is_same<ExprType1,ExprType2>::value, "" );

K_ * trans(K_);  // Results in a symmetric matrix
K_ * trans(L_);  // Does not result in a symmetric matrix

Whereas ExprType1 and ExprType2 are identical, the first expression will result in a symmetric matrix and the second will not. This shows that it would be necessary to determine symmetry at runtime, which unfortunately can result in a measurable overhead (especially for tiny matrices). Therefore it is always necessary to use declaim() manually.

I hope this helps to understand the rational behind declsym(),

Best regards,

Klaus!

[Mikhail Katliar]

Hi Klaus,

I think I understand the rationale behind declsym(), but my question is somewhat different. Since the SymmetricMatrix type restricts the set of possible values of instances by enforcing the a(i,j)=a(j,i) invariant, it does provide information about symmetry. Considering the following example,

using Matrix = blaze::DynamicMatrix<Real>;
using SymmetricMatrix = blaze::SymmetricMatrix<Matrix>;

SymmetricMatrix S;
Matrix K;
Matrix L;

K * trans(K);  // 1. Results in a symmetric matrix
K * trans(L);  // 2. Does not result in a symmetric matrix
S + S; // 3. Results in a symmetric matrix
K * S; // 4. Does not result in a symmetric matrix
S * trans(K); // 5. Does not result in a symmetric matrix
K * S * trans(K); // 6. Results in a symmetric matrix

I wonder why it is possible to detect symmetry at compile time in (1) and (3), but not in (6). Is it because (6) is a product of 3 expressions, and none of sub-products is symmetric?

[Klaus Iglberger]

Hi Mikhail!

I think there is a misunderstanding. Whereas Blaze can indeed determine that S+S will result in a symmetric matrix at compile time, Blaze cannot do that for any matrix multiplication (not even K_*trans(K_)). In my example I tried to point out that this is generally not possible because it depends on runtime information. The comment in the example should be adapted:

K_ * trans(K_);  // Results in a symmetric matrix at runtime, not at compile time; resulting type is 'DynamicMatrix<Real>'
declsym( K_ * trans(K_) );  // Results in a symmetric matrix at runtime and compile time; resulting type is 'SymmetricMatrix<DynamicMatrix<Real>>'

I apologize for my inaccurate reply,

Best regards,

Klaus!