How can I solve the sparse linear system with the coefficients in Z_2?

Question

I want to solve a large sparse linear equation systems with coefficients in Z_2 using Eigen. First we have tried the Boolean type which does not work because 1+1=1 in Boolean but I expect 1+1=0. Hence a solution might be a new customised scalar type. But how it works exactly? Some suggestions of other packages or software are also welcome.

Answer 1

As the operators for fundamental types can not be overloaded. You need a customized scalar type. Here's the basic document telling you how to do that.

http://eigen.tuxfamily.org/dox/TopicCustomizingEigen.html#CustomScalarType

Basically you need to do 3 things.

1. make sure the common operator (+,-,*,/,etc.) are supported by the type T
2. add a specialization of struct Eigen::NumTraits
3. define the math functions that makes sense for your type. This includes standard ones like sqrt, pow, sin, tan, conj, real, imag, etc, as well as abs2 which is Eigen specific. (see the file Eigen/src/Core/MathFunctions.h)

Actually you can only define those operators and math functions that are needed for equation solving.

The above link provides a simple example for type adtl::adouble. The code is very short because most of the operations are already well defined. In the Eigen source dir unsupported/, there's another 3rd-party type mpfr::mpreal support. You can start from this link.

https://eigen.tuxfamily.org/dox/unsupported/group__MPRealSupport__Module.html

In the Eigen source dir, these files are related to mpfr::mpreal support. They may be useful.

./unsupported/Eigen/MPRealSupport ./unsupported/test/mpreal/mpreal.h ./unsupported/test/mpreal_support.cpp

---

This is a minimum working example with matrix multiplication support in Z_2:

#include <iostream>
#include <Eigen/Eigen>

namespace Z2 {

struct Scalar {
  Scalar() :
      data(0) {
  }
  Scalar(bool x) :
      data(x) {
  }
  bool data;

  inline Scalar operator+=(const Scalar& a) {
    return data ^= a.data;
  }
};

inline Scalar abs(const Scalar& a) {
  return a;
}

inline Scalar operator+(const Scalar& a, const Scalar& b) {
  return a.data ^ b.data;
}

inline Scalar operator*(const Scalar& a, const Scalar& b) {
  return a.data & b.data;
}

template<class E, class CT>
std::basic_ostream<E, CT> &operator <<(std::basic_ostream<E, CT> &os,
                                       const Z2::Scalar &m) {
  return os << m.data;
}

}

namespace Eigen {
// follow all other traits of bool
template<> struct NumTraits<Z2::Scalar> : NumTraits<bool> {
  typedef Z2::Scalar Real;
  typedef typename internal::conditional<sizeof(Z2::Scalar) <= 2, float, double>::type NonInteger;
  typedef Z2::Scalar Nested;
};
}

int main(void) {
  using namespace Eigen;
  Matrix<Z2::Scalar, Dynamic, Dynamic> x(2, 2), y(2, 2), i(2, 2), j(2, 2);
  x.row(0) << 1, 1;
  y.col(0) << 1, 1;
  i.setIdentity();
  j = i.array() + 1;
  std::cout << "x=\n" << x << std::endl;
  std::cout << "y=\n" << y << std::endl;
  std::cout << "i=\n" << i << std::endl;
  std::cout << "j=\n" << j << std::endl;
  std::cout << "x+y=\n" << x + y << std::endl;
  std::cout << "x.*y=\n" << x.array() * y.array() << std::endl;
  std::cout << "y*j=\n" << y * j << std::endl;
  std::cout << "abs(x)=\n" << x.array().abs() << std::endl;
  return 0;
}

Result:

x=
1 1
0 0
y=
1 0
1 0
i=
1 0
0 1
j=
0 1
1 0
x+y=
0 1
1 0
x.*y=
1 0
0 0
y*j=
0 1
0 1
abs(x)=
1 1
0 0