Hermitian-ness is not preserved when adding identity matrix

[pkerichang]

Simple code to reproduce the bug:

#include <iostream>

#include <blaze/math/DynamicMatrix.h>
#include <blaze/math/HermitianMatrix.h>
#include <blaze/math/IdentityMatrix.h>

int main(int argc, char **argv) {

  using mat_t =
      blaze::HermitianMatrix<blaze::DynamicMatrix<blaze::complex<double>>>;
  using iden_t = blaze::IdentityMatrix<blaze::complex<double>>;

  int N = 4;
  mat_t A;
  resize(A, N, N);
  randomize(A);

  iden_t I1(N);

  auto tmp1 = 3.5 * A - 4.1 * I1;

  std::cout << "tmp1 expr is hermitian: "
            << blaze::IsHermitian_v<decltype(tmp1)> << std::endl;

  mat_t I2(I1);

  auto tmp2 = 3.5 * A - 4.1 * I2;

  std::cout << "tmp2 expr is hermitian: "
            << blaze::IsHermitian_v<decltype(tmp2)> << std::endl;

  return 0;
}

terminal output:

tmp1 expr is hermitian: 0
tmp2 expr is hermitian: 1

Seems like for some reason the code did not think IdentityMatrix is hermitian.

One thing to note is that this bug rarely impacts behavior, because when an expression is assigned to a hermitian matrix, if IsHermitian_v<> returns False, the matrix will actually check the content of the matrix and if it’s hermitian, the expression will get assigned successfully. However, this bug pop up for me because I’m working with matrices with large variations in magnitudes, and apparently the numeric noise is large enough for the content check to fail. This brings up an orthogonal question of mine: is there a way to set the relative/absolute error tolerances of floating point equality checks in blaze?

My current workaround is to save the identity matrix inside a hermitian matrix. However, this may waste unncessary memory.

[Klaus Iglberger]

Hi!

Thanks a lot for taking the time to report this potential issue. You are correct, the expression 4.1 * I1 does not preserve the Hermitian property of the identity matrix. This is not a bug, though, but a limitation of the currently available types in Blaze: There is no type, that would be able to express that the result is both diagonal and Hermitian, but not symmetric or identity.

For that reason Blaze provides declaration operations. In order to declare an expression as Hermitian you can use the declherm() expression:

auto tmp1 = declherm( 3.5 * A - 4.1 * I1 );
std::cout << "tmp1 expr is hermitian: "
            << blaze::IsHermitian_v<decltype(tmp1)> << std::endl;

Output:

tmp1 expr is hermitian: 1

I hope this explains the rational and helps to deal with this situation.

Best regards,

Klaus!

[pkerichang]

Hello,

Understood. However, this difficulty came up when I was working on pull request #49. Consider the following code example:

template <typename Mat>
void a_function(const Mat A) {
    IdentityMatrix<ElementType_t<Mat>> I( A.rows() );
    auto tmp = 3.5 * A + 4.1 * I;
    // do more things with tmp
    // ...
}

here, if A is upper/lower triangular, symmetric, or hermitian, then we know tmp will also keep those adaptor types (whereas for strictly lower/upper triangular adaptors, tmp will lose them). However, we cannot use any of the declX() functions because the adaptor type of A (if any) is not known.

if we allow C++17 then this problem can be solved with constexpr-if, although the solution will involve bunch of if statements and looks ugly. Alternatively, Is it possible to write a templated generic decl() function with signature like:

template <typename M, 
          bool isUpper, 
          bool isLower, 
          bool isSymmetric, 
          bool isHermitian, 
          bool isStrictUpper, 
          bool isStrictLower> 
auto decl(M x);

that basically turns into one of the declX() functions based on the boolean flags. Then you can write something like:

auto tmp = decl<IsUpper_v<M>, IsLower_v<M>, IsSymmetric_v<M>, IsHermitian_v<M>, false, false>( 3.5 * A + 4.1 * I );

to achieve our goal.

If the multiple booleans is undesirable (as it is hard to extend new adaptor types in the future), we can assign a unique integer to every adaptor type, then you can have functions and type traits structs like:

template <int AdaptorID, typename M> 
auto decl(M x);

template <class M> struct AdaptorID;
template <class M> inline constexpr int AdaptorID_v = AdaptorID<M>::value;  

so you write code like:

auto tmp = decl<AdaptorID_v<Mat>>( 3.5 * A + 4.1 * I );

to achieve the functionality of dropping strictly lower/upper triangular types, you can achieve it with adding helper functions that manipulate the adaptor IDs, so the code becomes:

constexpr int dropStrictLower(int id);
constexpr int dropStrictUpper(int id);

auto tmp = decl<dropStrictLower(dropStrictUpper(AdaptorID_v<Mat>))>( 3.5 * A + 4.1 * I );

‌

In summary, I think if there’s a templated generic decl method, this will help a lot in writing generate templated code that is meant to work with multiple adaptor types.

Do you think this is a feature that sounds interesting?

[Klaus Iglberger]

Hi!

Thanks for the additional details. I completely understand the problem. Allow me to investigate this in the context of pull request #49.

Best regards,

Klaus!