Gauss-Bareissによる整数行列の行列式

README.md

Gauss-Bareissによる整数行列の行列式

整数行列の行列式を求めるための素朴な方法としては有理数行列のLU分解が考えられるが、これは計算途中で分子分母が大きくなるため効率が悪い。

整数行列の行列式を求める効率的な方法として以下の2つが考えられる:

• 十分多くのpについてLU分解等でZ/pZ上の行列式を計算し、Chinese Remainder Theoremで復元する。
• 整数計算のみで掃き出しを行い、ちょうど割り切れることが保証されている除算によって計算を効率化する (Gauss-Bareiss)

ここでは後者を実装してみた。

例題の説明

例題は、0,1,-1からなる歪対称行列の行列式の平方根である。整数値歪対称行列の行列式の平方根は整数になり、Pfaffianと呼ばれる。Pfaffianは例えば平面グラフの完全マッチングの個数(FKT Algorithm)で使われる。mod nでの平方根は一意に決まらないから、mod nでの計算によってではこの値は求まらない。したがってこの問題は整数行列の行列式のよい例題となっている。

実装

• 言語 : C++11

• 行列 : Eigen3

• 多倍長整数 : gmpxx

• gb.cpp : Gauss-Bareissによる手法

• lu.cpp : 素朴な方法(有理数でのLU分解？)

• lu-double.cpp : 倍精度浮動小数点数での近似

実行方法 : 第1引数(任意)が行列の大きさ、第2引数(任意)が乱数のシード

実行結果

n=50

• gb : 0.02s
• lu : 0.04s

n=100

• gb : 0.11s
• lu : 0.46s

n=150

• gb : 0.31s
• lu : 2.27s

n=200

• gb : 0.92s
• lu : 7.79s

n=250

• gb : 2.35s
• lu : 21.64s

n=300

• gb : 4.59s
• lu : 47.53s

gb.cpp

#include <iostream>
#include <random>
#include <string>
#include <sstream>
#include <Eigen/Dense>
#include <gmpxx.h>

template<typename Scalar, int Rows, int Cols, int Options,
         int MaxRows, int MaxCols>
Scalar gb_determinant(
  Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> a) {
  int n = a.rows();
  assert(n == a.cols());
  int sign = 1;
  Scalar c = 1;
  for(int k = n-1; k >= 0; --k) {
    Scalar p = a(k, k);
    if(p == 0) {
      int k2;
      for(k2 = k-1; ; --k2) {
        if(k2 < 0) return 0;
        if(a(k2, k) != 0) break;
      }
      a.row(k2).swap(a.row(k));
      p = a(k, k);
      sign = -sign;
    }
    for(int i = 0; i < k; ++i) {
      for(int j = 0; j < k; ++j) {
        a(i, j) = (p * a(i, j) - a(i, k) * a(k, j)) / c;
      }
    }
    a.conservativeResize(k, k);
    c = p;
  }
  return sign * c;
}

int main(int argc, char *argv[]) {
  int n = 10;
  if(argc > 1) {
    std::istringstream is(argv[1]);
    is >> n;
  }
  int rand_seed = 12345;
  if(argc > 2) {
    std::istringstream is(argv[2]);
    is >> rand_seed;
  }
  std::mt19937 rand_source(rand_seed);
  auto binrand =
    std::bind(std::uniform_int_distribution<int>(0, 1),
              std::ref(rand_source));
  Eigen::Matrix<mpz_class, Eigen::Dynamic, Eigen::Dynamic> a(n, n);
  for(int i = 0; i < n; ++i) {
    for(int j = 0; j < n; ++j) {
      a(i, j) = binrand();
    }
  }
  Eigen::Matrix<mpz_class, Eigen::Dynamic, Eigen::Dynamic> b
    = a - a.transpose();
  mpz_class det_b = gb_determinant(b);
  mpz_class pfaffian = sqrt(det_b);
  std::cout << pfaffian << std::endl;
  return 0;
}

lu-double.cpp

#include <iostream>
#include <random>
#include <string>
#include <sstream>
#include <Eigen/Dense>

int main(int argc, char *argv[]) {
  int n = 10;
  if(argc > 1) {
    std::istringstream is(argv[1]);
    is >> n;
  }
  int rand_seed = 12345;
  if(argc > 2) {
    std::istringstream is(argv[2]);
    is >> rand_seed;
  }
  std::mt19937 rand_source(rand_seed);
  auto binrand =
    std::bind(std::uniform_int_distribution<int>(0, 1),
              std::ref(rand_source));
  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> a(n, n);
  for(int i = 0; i < n; ++i) {
    for(int j = 0; j < n; ++j) {
      a(i, j) = binrand();
    }
  }
  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> b
    = a - a.transpose();
  double det_b = b.determinant();
  double pfaffian = sqrt(det_b);
  std::cout << pfaffian << std::endl;
  return 0;
}

lu.cpp

#include <iostream>
#include <random>
#include <string>
#include <sstream>
#include <Eigen/Dense>
#include <gmpxx.h>

int main(int argc, char *argv[]) {
  int n = 10;
  if(argc > 1) {
    std::istringstream is(argv[1]);
    is >> n;
  }
  int rand_seed = 12345;
  if(argc > 2) {
    std::istringstream is(argv[2]);
    is >> rand_seed;
  }
  std::mt19937 rand_source(rand_seed);
  auto binrand =
    std::bind(std::uniform_int_distribution<int>(0, 1),
              std::ref(rand_source));
  Eigen::Matrix<mpq_class, Eigen::Dynamic, Eigen::Dynamic> a(n, n);
  for(int i = 0; i < n; ++i) {
    for(int j = 0; j < n; ++j) {
      a(i, j) = binrand();
    }
  }
  Eigen::Matrix<mpq_class, Eigen::Dynamic, Eigen::Dynamic> b
    = a - a.transpose();
  mpz_class det_b = b.determinant().get_num();
  mpz_class pfaffian = sqrt(det_b);
  std::cout << pfaffian << std::endl;
  return 0;
}

Makefile

#!/usr/bin/make -f

CXX = g++
CXXFLAGS = -std=c++11 -O2 -Wall -Wextra -g
CPPFLAGS = -I /usr/include/eigen3/
LDLIBS = -lm -lgmpxx -lgmp

EXECS = gb lu lu-double

all: $(EXECS)

clean:
	$(RM) $(EXECS)