定義に従って行列式を求める

ここでいう行列式の定義とはこれです。 \[{\rm det}A=\sum_{\sigma \in S_n} {\rm sgn}(\sigma)\prod_{i=1}^n a_{i \sigma(i)}\] ここで言う$\sigma$は置換で、すべての置換の総和で行列式を定義したものですね。ちなみに${\rm sgn}$は偶置換、奇置換それぞれで1,-1を取る関数です。詳しくは線形代数の教科書でも見てください。

この定義、何が恐ろしいって計算量が${\rm O}(n \times n!)$なんですよね。すべての置換パターンで積を計算しなければならないので、n個の数の順列分の計算量が生まれるというわけです。行列の次数が上がるにつれて指数オーダーよりも激しく計算量が増えていきます。
ただ、それぞれの置換パターンは独立していますので並列処理も可能です。ですので、ベンチマークに使いやすいわけです。

主たる課題は下記のとおりです。

• n次の行列のときnつの数の積を取るためオーバーフローしやすい（＝多倍長整数型などを作る必要がある）
• すべての置換パターンを求める必要がある
• 並列処理の効率化（クリティカルセクション等々）

大学1年生のときC言語で計算するプログラムを書いて遊んでいましたが、なかなか大変でした。でも、C#は様々なライブラリがあり、並列処理も簡単に書けます。というわけで、C#で書いてみました。
コンピューターでの小数には「精度」が伴ってしまいますので、行列の要素は整数に限定してプログラムを作っていきます。

多倍長整数

C#にはintやlongなどの型があります。それぞれ32bit、64bitとビット数は決まっており、それより巨大な計算結果になるような演算をすると正確な値が得られなくなってしまいます。

そこで、必要に応じてビット数を拡張し、特にビット数の縛りなくメモリが許す限りどんな大きな整数でも扱えるような機構が欲しくなります。C#のような言語ならばそういったクラスや構造体を作れば良でしょう。
というのは誰もが思うことで、すでにそういうのが用意されています。

BigInteger 構造体

使うときに参照を追加してやる必要がありますが、それだけで簡単に多倍長整数で計算してくれます。演算子もオーバーロードされており、特に違和感なく計算ができるところがとても良いです。

ですが、さすがにLINQのメソッドまでは用意されておりませんので、今回のプログラムで使う分を実装しておきましょう。

public static class BigIntegerEnumerable
{
    public static IEnumerable<BigInteger> Range(BigInteger start, BigInteger count)
    {
        for(BigInteger i = 0; i < count; i++)
            yield return start + i;
    }
}

public static class BigIntegerUtil
{
    public static BigInteger Factorial(int start) => Enumerable.Range(1, start).Aggregate(new BigInteger(1), (total, now) => total * now);

    public static BigInteger Sum(this IEnumerable<BigInteger> source) => source.Aggregate(BigInteger.Zero, (now, next) => now + next);

    public static BigInteger Product(this IEnumerable<BigInteger> source) => source.Aggregate(BigInteger.One, (now, next) => now * next);
}

置換は(1, 2, ..., n)と並んでいる数列を入れ替えていくという話なので、Enumerable.Rangeみたいなメソッドがあるととても便利です。
また、多倍長整数を使うモチベーションとして「階乗」というのは非常に大きいはずですが、なぜかBigInteger構造体には階乗を求めるメソッドが無いので実装しておきます。あとは総和と総乗を取るのでSumとProductという拡張メソッドを作っておきます。Aggregate大活躍です。

順列と置換

例えば(1,2,3)の3つの数の順列は、(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1)の6パターンあります。そして、順列を小さい順？に並べるとこの順序になるはずです。
じゃあ任意の数列に対してi番目の並び替えは何か、というのを求めるメソッドを用意します。

static int[] Permutation(int[] source, BigInteger index)
{
    if(source.Length <= 1)
        return source.ToArray();

    var list = source.ToList();
    var div = (int)BigInteger.DivRem(index, BigIntegerUtil.Factorial(source.Length - 1), out BigInteger remainder);

    var first = list[div];
    list.RemoveAt(div);

    return new int[] { first }.Concat(Permutation(list.ToArray(), remainder)).ToArray();
}

今の配列の一番前に来るべき数を決めて、残りの項目を1つの数列として再帰的にこの関数を呼び出すことでindex番目の並び順を求めています。

ただ、このindexと置換の回数は直接的な相関はありません（たぶん）。ですので、置換の回数を求めるメソッドも用意しました。

static int ReplaceCount(int[] array1, int[] array2)
{
    if(array1.Length != array2.Length)
        throw new ArgumentException("Invalid array lengthes");

    int count = 0;            

    for(int i = 0; i < array1.Length; i++) {
        if(array1[i] != array2[i]) {
            var nextIndex = array2.ToList().IndexOf(array1[i], i + 1);
            if(nextIndex < 0)
                throw new InvalidOperationException("Cannot find same value.");

            var swap = array2[i];
            array2[i] = array2[nextIndex];
            array2[nextIndex] = swap;
            count++;
        }
    }

    return count;
}

やってることは選択ソートみたいなもんです。選択ソートについて置換の回数を見ているので、別のアルゴリズムを使えば回数は異なる可能性があります。ただ、偶奇性は一意に定まります。

並列実装

ここまで来れば、定義に従って行列式を求めるだけです。

Random rnd = new Random();
int[,] matrix = new int[order, order];

for(int i = 0; i < order; i++) {
    for(int j = 0; j < order; j++)
        matrix[i, j] = rnd.Next(10);
}

for(int j = 0; j < order; j++) {
    for(int i = 0; i < order; i++)
        Console.Write($"{matrix[i, j]} ");
    Console.WriteLine();
}

Dictionary<int, BigInteger> threaddet = new Dictionary<int, BigInteger>();
object threadlock = new object();

Stopwatch sw = new Stopwatch();
sw.Start();

var source = Enumerable.Range(0, order).ToArray();
Parallel.ForEach(BigIntegerEnumerable.Range(0, BigIntegerUtil.Factorial(order)), patternIndex => {
    var pattern = Permutation(source, patternIndex);
    var res = source.Zip(pattern, (r, c) => (BigInteger)matrix[r, c]).Product() * (ReplaceCount(source, pattern) % 2 == 0 ? 1 : -1);

    var id = System.Threading.Thread.CurrentThread.ManagedThreadId;
    if(!threaddet.ContainsKey(id)) {
        lock(threadlock)
            threaddet.Add(id, 0);
    }
    threaddet[id] += res;
});

var det = threaddet.Values.Sum();

sw.Stop();

Console.WriteLine($"det = {det}");
Console.WriteLine($"Elapsed: {sw.Elapsed.TotalSeconds}s");

まずはorder次の正方行列をランダムに生成します。わかりやすく各要素は0～9の整数にしています。
元の行列が決まったらorderの階乗だけ置換パターンがあるのでそれをParallel.ForEachで並列処理をしています。それぞれのパターンについて要素を取ってきて積を取ってから符号をつけます。
積が求まったら総和を取らなければなりませんが、和を取るためにはスレッド間での同期をとらなければなりません。ですので、それは一番最後に回すことにしています。具体的にはそれぞれのスレッドで計算した分の総和のみを求めており、すべてのスレッドの処理が終わったら各スレッドの総和を足し合わせる、という処理になっています。

スレッドの総和を求めるためには、スレッドIDをキーとしたDictionaryを用意して求めます。キーの追加時はロックをしないと、おそらく内部的なリストの「最後尾」にキーを追加というタイミングで複数スレッド間での競合が起きバグります。

ちゃんと検算もしたで pic.twitter.com/EOsAuwrok6

— EH500金太郎 (@EH500_Kintarou) 2017年11月18日

私のPCはi7-4790Kを搭載していますが、11次の行列で2分半ちょっと計算にかかりました。行列を計算し始めてからOctaveを起動して行列を打ち込んで行列式を求めてもまだ定義に従った計算では求まらない、といったくらいのペースになります。

昔、C言語で作った時はえらく苦労をした記憶がありますが、C#だとこんなにすんなりできてしまいました。
やはり、並列化が体系化されていたり、演算子オーバーロード等で非常に使いやすい多倍長整数型が用意されていたり、あとはLINQで配列を簡単に扱えたりするのはとても大きいですね。

もうちょっとアルゴリズムが改良できるぞ、とかあったら是非教えてください。
特に置換パターンを求めたり置換回数を求めるあたりはまだ何かやりようがある気がしています。あ、「掃き出し法を使うと速いぞ」とか言うのは無しで。あくまでも定義に従った計算方式ということで。