第22章
行列式を計算する三つの方法
（やさしい版）
Pearls of Functional Algorithm Design（関数プログラミングによるアルゴリズム設計の真珠）
どんな問題？
n × n の行列 A の行列式 |A| を、コンピュータでうまく計算したい、というのが今回のテーマです。
行列式の定義自体は「すべての置換について符号付きの積を足す」というライプニッツの公式で与えられ、次のように書けます。
|A| = Σπ sign (π) Π1 ≤ j ≤ n ajπ(j)
ただし、この定義をそのまま実装すると、n! に比例する時間がかかってしまいます（n = 10 でも約 360 万通り、n = 20 なら宇宙の年齢でも足りません）。もっと現実的な方法が必要です。
この章で扱う三つの方法
- 方法1：ガウスの消去法（有理数除算を使う）
- 方法2：キオの枢軸凝縮法（整数除算だけで済ませる）
- 方法3：行列積の反復（除算をまったく使わない）
それぞれ得意・不得意があり、最後に実測で比べます。
教科書的な方法（ウォームアップ）
小行列式を再帰で展開する
まず、高校で習う余因子展開を素直に書いてみましょう。例えば 3 × 3 の場合はこうです。
| a11 a12 a13 |
| a21 a22 a23 | = a11 | a22 a23 | − a21 | a12 a13 | + a31 | a12 a13 |
| a31 a32 a33 | | a32 a33 | | a32 a33 | | a22 a23 |
行列を「行のリスト」で表すと、Haskell では次のように書けます。
det :: [[Integer]] → Integer
det [[x]] = x
det xss = foldr1 (−) (zipWith (∗) col1 (map det (minors cols)))
where col1 = map head xss
cols = map tail xss
Dart
// 余因子展開による行列式（Θ(n!)）
BigInt detCofactor(List<List<BigInt>> xss) {
if (xss.length == 1 && xss[0].length == 1) return xss[0][0];
final col1 = xss.map((row) => row.first).toList(); // 第 1 列
final cols = xss.map((row) => row.sublist(1)).toList(); // 第 1 列を除いた行列
final subDets = minors(cols).map(detCofactor).toList(); // 各小行列式
// foldr1 (-) で交互に加減算： a0 - (a1 - (a2 - ...))
BigInt acc = col1.last * subDets.last;
for (int i = col1.length - 2; i >= 0; i--) {
acc = col1[i] * subDets[i] - acc;
}
return acc;
}
1 × 1 の場合はその値そのもの。それ以外は「第 1 列の各要素 × 対応する小行列の行列式」を、符号を交互にしながら足し引きします。
minors は「先頭の要素を 1 つずつ抜いた列（小行列に対応）」を作る関数です。
minors :: [a] → [[a]]
minors [] = []
minors (x : xs) = xs : map (x :) (minors xs)
Dart
// minors [a,b,c,d] = [[b,c,d], [a,c,d], [a,b,d], [a,b,c]]
List<List<T>> minors<T>(List<T> xs) {
if (xs.isEmpty) return <List<T>>[];
final x = xs.first;
final rest = xs.sublist(1);
final tail = minors(rest); // 再帰
return [rest, ...tail.map((ys) => [x, ...ys])];
}
例：minors "abcd" = ["bcd", "acd", "abd", "abc"]。
問題点
短くて可愛いコードですが、計算量は Θ(n!)。n が 2 か 3 くらいなら十分ですが、それ以上になるとまったく実用になりません。だからこれから紹介する 3 つの方法が必要になるわけです。
方法1：有理数除算を使う（ガウスの消去法）
基本アイデア
ガウスの消去法は、次の性質を使います。
- ある行の何倍かを別の行に足しても、行列式の値は変わらない。
これを使って、第 1 行の適当な倍数を他の行に足して第 1 列を全部 0 にします。次は残った (n−1) × (n−1) の部分行列に対して同じことを繰り返し、最終的に上三角行列に変形します。上三角行列の行列式は対角要素の積です。
ピボットに注意
ややこしいのは、先頭要素が 0 になっているケース。この場合は、先頭要素が 0 でない行（これをピボットと呼びます）を探して、それを第 1 行と入れ替える必要があります。行の入れ替えは、位置が奇数回動いた場合に行列式の符号を反転させます。
det :: [[Ratio Integer]] → Ratio Integer
det [[x]] = x
det xss =
case break ((≠ 0) · head) xss of
(yss, []) → 0
(yss, zs : zss) → let x = head zs ∗ det (reduce zs (yss ++ zss))
in if even (length yss) then x else − x
Dart
// 方法1: ガウスの消去法（有理数）。Rational は Dart には無いので簡易分数型を用意。
class Rat {
final BigInt n, d; // 分子・分母（既約）
Rat._(this.n, this.d);
factory Rat(BigInt n, [BigInt? d]) {
d ??= BigInt.one;
if (d == BigInt.zero) throw ArgumentError('zero denom');
final g = n.gcd(d).abs();
var nn = n ~/ g, dd = d ~/ g;
if (dd.isNegative) { nn = -nn; dd = -dd; }
return Rat._(nn, dd);
}
static final zero = Rat(BigInt.zero);
Rat operator +(Rat o) => Rat(n * o.d + o.n * d, d * o.d);
Rat operator -(Rat o) => Rat(n * o.d - o.n * d, d * o.d);
Rat operator *(Rat o) => Rat(n * o.n, d * o.d);
Rat operator /(Rat o) => Rat(n * o.d, d * o.n);
Rat operator -() => Rat(-n, d);
bool get isZero => n == BigInt.zero;
}
Rat detGauss(List<List<Rat>> xss) {
if (xss.length == 1 && xss[0].length == 1) return xss[0][0];
// 先頭要素が 0 でない行を探す（break ((≠ 0) · head) 相当）
final k = xss.indexWhere((row) => !row.first.isZero);
if (k < 0) return Rat.zero;
final zs = xss[k];
final rest = [...xss.sublist(0, k), ...xss.sublist(k + 1)];
final x = zs.first * detGauss(reduce(zs, rest));
return k.isEven ? x : -x; // 行入れ替えの符号
}
break でリストを「先頭が 0 の部分（yss）」と「先頭が 0 でない部分（zs : zss）」に分け、後者が空なら行列は特異（行列式は 0）。そうでなければ、zs をピボットとして残りの行を簡約します。
reduce xs yss = map (reduce1 xs) yss
reduce1 (x : xs) (y : ys) = zipWith (λa b → b − d ∗ a) xs ys
where d = y / x
Dart
// ピボット行 xs で他の行たち yss の第 1 列を 0 にする（先頭も落として (n-1)x(n-1) に）
List<List<Rat>> reduce(List<Rat> zs, List<List<Rat>> yss) =>
yss.map((y) => reduce1(zs, y)).toList();
List<Rat> reduce1(List<Rat> xs, List<Rat> ys) {
final d = ys.first / xs.first; // 消去係数
return [
for (var i = 1; i < xs.length; i++) ys[i] - d * xs[i],
];
}
短所
除算に有理数（分数）を使うため、毎回「分子と分母の最大公約数を計算して約分」しなくてはならず、そこに時間を取られます。整数のまま扱えると嬉しいのですが……。
方法2：整数除算を使う（キオの枢軸凝縮法）
キオの恒等式
a11 ≠ 0 のとき、次の等式が成り立ちます（これがキオの恒等式）。
|A| = |X| / a11n−2
ここで X は (n−1) × (n−1) の「凝縮した行列」で、各要素は次の 2 × 2 行列式で決まります。
xjk = | a11 a1k |
| aj1 ajk |
ポイント
この方法でも除算は必要ですが、必ず割り切れることが保証されているので、有理数ではなく整数除算で済みます。約分の手間がありません。
実装
a11 が 0 なら、ガウスの消去法と同じくピボットを探して行を交換します。
det :: [[Integer]] → Integer
det [[x]] = x
det xss =
case break ((≠ 0) · head) xss of
(yss, []) → 0
(yss, zs : zss) → let x = det (condense (zs : yss ++ zss))
d = head zs ↑ (length xss − 2)
y = x div d
in if even (length yss) then y else − y
Dart
// 方法2: キオの枢軸凝縮法（整数除算のみ）
BigInt detChio(List<List<BigInt>> xss) {
if (xss.length == 1 && xss[0].length == 1) return xss[0][0];
final k = xss.indexWhere((row) => row.first != BigInt.zero);
if (k < 0) return BigInt.zero;
final zs = xss[k];
// ピボット行を先頭に持ってきた行列を渡して凝縮
final swapped = [zs, ...xss.sublist(0, k), ...xss.sublist(k + 1)];
final x = detChio(condense(swapped));
final d = zs.first.pow(xss.length - 2); // a11^(n-2)
final y = x ~/ d; // 割り切れる
return k.isEven ? y : -y;
}
(↑) は指数演算です。condense は「第 1 行と他の各行をペアにして、それらの 2 × 2 行列式を並べる」処理を行います。
condense = map (map det · pair · uncurry zip) · pair
where pair (x : xs) = map ((,) x) xs
det ((a, b), (c, d)) = a ∗ d − b ∗ c
Dart
// 第 1 行 zs と他の各行 r をペアにして、列同士も第 1 列と組み、2x2 行列式を並べる
List<List<BigInt>> condense(List<List<BigInt>> xss) {
final zs = xss.first; // 第 1 行
final rows = xss.sublist(1); // 残りの行
return [
for (final r in rows) [
for (var k = 1; k < zs.length; k++)
zs[0] * r[k] - zs[k] * r[0], // 2x2 行列式
],
];
}
短所
凝縮 1 回で Θ(n2) 手間、それを n 回繰り返すので全体は Θ(n3)。有理数は使いませんが、代わりに途中の整数がどんどん巨大化するのが困りものです。
方法2の改良：凝縮と除算を交互に
交互にやれば数字が小さいまま
実は行列式には次の性質があります。
- A から得た凝縮行列 X、さらにその凝縮行列 Y について、Y の各要素は必ず a11 で割り切れる（a11 ≠ 0 のとき）。
この性質を使うと、キオの方法の最後にまとめて割っていた a11n−2 を、凝縮のたびに少しずつ割る形に変えられます。すると det = det′ 1 として、以下のように書けます。
det′ :: Integer → [[Integer]] → Integer
det′ k [[x]] = x
det′ k xss =
case break ((≠ 0) · head) xss of
(yss, []) → 0
(yss, zs : zss) → let x = det′ (head zs) (cd k (zs : yss ++ zss))
in if even (length yss) then x else − x
Dart
// 改良版キオ: 凝縮と除算を交互に。det(A) = detPrime(1, A)
BigInt detPrime(BigInt k, List<List<BigInt>> xss) {
if (xss.length == 1 && xss[0].length == 1) return xss[0][0];
final j = xss.indexWhere((row) => row.first != BigInt.zero);
if (j < 0) return BigInt.zero;
final zs = xss[j];
final swapped = [zs, ...xss.sublist(0, j), ...xss.sublist(j + 1)];
final x = detPrime(zs.first, cd(k, swapped)); // 次の除数は今の a11
return j.isEven ? x : -x;
}
BigInt detChioImproved(List<List<BigInt>> xss) => detPrime(BigInt.one, xss);
cd（condense and divide の略）は「凝縮したら同時に k で割る」処理です。
cd k = map (map det · pair · uncurry zip) · pair
where pair (x : xs) = map ((,) x) xs
det ((a, b), (c, d)) = (a ∗ d − b ∗ c) div k
Dart
// condense と同じだが、各要素を k で割る（割り切れる）
List<List<BigInt>> cd(BigInt k, List<List<BigInt>> xss) {
final zs = xss.first;
final rows = xss.sublist(1);
return [
for (final r in rows) [
for (var j = 1; j < zs.length; j++)
(zs[0] * r[j] - zs[j] * r[0]) ~/ k,
],
];
}
もちろん、除算の回数自体は Θ(n3) まで増えます。でも整数が小さく保たれる分、全体としてはずっと有利です。
方法3：除算を使わない（反復行列積）
MUT という魔術的な行列
最後に、除算をまったく使わない方法を紹介します。少し不思議なやり方で、証明はここでは省略します。
n × n 行列 X に対して、MUT(X)（make upper triangular の略）を次のように定義します。
- 対角より下の要素はすべて 0 にする。
- 対角より上の要素はそのまま残す。
- 対角要素 i 番目は、その位置より下側にある対角要素の総和の符号を反転したもので置き換える。
MUT(X) =
( −Σj=2n xjj x12 ... x1n )
( 0 −Σj=3n xjj ... x2n )
( ... )
( 0 0 ... −Σj=n+1n xjj )
（最後の対角要素は、それより下がないので 0 になります。）
反復で |A| が出てくる
FA(X) = MUT(X) × A と定義し、B = FAn−1 A′ と置きます（n が奇数なら A′ = A、偶数なら A′ = −A）。すると、この B はほとんど 0 の行列で、左上の b11 だけが |A| に等しくなります。
計算量
MUT × A の計算に Θ(n3)、それを n − 1 回繰り返すので、全体で Θ(n4) です。除算は要らないけれど、乗算は増えます。
実装
det :: [[Integer]] → Integer
det ass = head (head bss)
where
bss = foldl (matmult · mut) ass′ (replicate (n − 1) ass)
ass′ = if odd n then ass else map (map negate) ass
n = length ass
Dart
// 方法3: 反復行列積による行列式（除算なし）
BigInt detMatmult(List<List<BigInt>> ass) {
final n = ass.length;
// n が偶数なら符号反転から始める
var bss = n.isOdd
? ass
: ass.map((r) => r.map((e) => -e).toList()).toList();
// n-1 回、bss ← matmult(mut(bss), ass)
for (var i = 0; i < n - 1; i++) {
bss = matmult(mut(bss), ass);
}
return bss.first.first;
}
mut は MUT の実装で、対角の下側に 0 を敷き詰め、対角要素をずらした累積和で置き換えます。
mut xss
= zipWith (++) zeros (zipWith (:) ys (zipWith drop [1..] xss))
where ys = map negate (tail (scanr (+) 0 (diagonal xss)))
Dart
// MUT: 対角線より下は 0、対角要素は「その下の対角要素の総和」の符号反転、上はそのまま
List<List<BigInt>> mut(List<List<BigInt>> xss) {
final n = xss.length;
final diag = diagonal(xss);
// scanr (+) 0 diag = [d0+d1+...+dn-1, d1+...+dn-1, ..., dn-1, 0]
// ys = -tail(scanr) = [-(d1+...+dn-1), -(d2+...+dn-1), ..., -dn-1, 0]
final ys = <BigInt>[];
var s = BigInt.zero;
for (var i = n - 1; i >= 1; i--) {
s += diag[i];
ys.insert(0, -s);
}
ys.add(BigInt.zero); // 最後の対角要素
return [
for (var i = 0; i < n; i++) [
for (var j = 0; j < i; j++) BigInt.zero, // 対角より下
ys[i], // 対角
...xss[i].sublist(i + 1), // 対角より上
],
];
}
zeros = [take j (repeat 0) | j ← [0..]]
Dart
// zeros = [[], [0], [0,0], [0,0,0], ...]（無限リスト）
Iterable<List<BigInt>> zeros() sync* {
for (var j = 0; ; j++) {
yield List<BigInt>.filled(j, BigInt.zero);
}
}
diagonal [] = []
diagonal (xs : xss) = head xs : diagonal (map tail xss)
Dart
// 対角成分の抽出： xss[i][i] を並べる
List<BigInt> diagonal(List<List<BigInt>> xss) =>
[for (var i = 0; i < xss.length; i++) xss[i][i]];
matmult xss yss = zipWith (map · dp) xss (repeat (transpose yss))
dp xs ys = sum (zipWith (∗) xs ys)
Dart
// 内積
BigInt dp(List<BigInt> xs, List<BigInt> ys) {
var s = BigInt.zero;
for (var i = 0; i < xs.length; i++) s += xs[i] * ys[i];
return s;
}
// 行列積
List<List<BigInt>> matmult(List<List<BigInt>> xss, List<List<BigInt>> yss) {
final cols = transpose(yss);
return [for (final row in xss) [for (final col in cols) dp(row, col)]];
}
List<List<T>> transpose<T>(List<List<T>> m) {
if (m.isEmpty) return <List<T>>[];
final w = m[0].length;
return [for (var j = 0; j < w; j++) [for (final row in m) row[j]]];
}
上三角性を活かした最適化
MUT(X) は「対角より下の要素」に依存しません。遅延評価のおかげで実際にはそれらは計算されませんが、それでも「上三角行列と任意の行列の積で、結果も上三角行列にする」専用の乗算 trimult を作る方が効率的です。
trimult xss yss = zipWith (map · dp) xss (submats (transpose yss))
Dart
// 上三角 xss と yss の積で、結果も上三角にする。i 行目は yss^T の第 i 部分行列を使う
List<List<BigInt>> trimult(List<List<BigInt>> xss, List<List<BigInt>> yss) {
final subs = submats(transpose(yss)); // 上から順に取る部分行列
return [
for (var i = 0; i < xss.length; i++)
[for (final col in subs[i]) dp(xss[i].sublist(i), col)],
];
}
submats :: [[a]] → [[[a]]]
submats [[x]] = [[[x]]]
submats xss = xss : submats (map tail (tail xss))
Dart
// 主部分行列の列： M, M の左上を 1 行 1 列削ったもの, さらに削ったもの, ...
List<List<List<T>>> submats<T>(List<List<T>> xss) {
final result = <List<List<T>>>[];
var m = xss;
while (true) {
result.add(m);
if (m.length == 1 && m[0].length == 1) break;
m = [for (var i = 1; i < m.length; i++) m[i].sublist(1)];
}
return result;
}
上三角行列 xss に対しては、mut の定義もぐっとシンプルになります。
mut xss = zipWith (:) ys (map tail xss)
where ys = map negate (tail (scanr (+) 0 (map head xss)))
Dart
// 上三角行列専用の MUT。各行はすでに「対角以降」だけを保持している前提。
// map head xss がその行の対角要素になる。
List<List<BigInt>> mutTri(List<List<BigInt>> xss) {
final heads = [for (final row in xss) row.first]; // 対角要素
final ys = <BigInt>[];
var s = BigInt.zero;
for (var i = heads.length - 1; i >= 1; i--) {
s += heads[i];
ys.insert(0, -s);
}
ys.add(BigInt.zero);
return [for (var i = 0; i < xss.length; i++) [ys[i], ...xss[i].sublist(1)]];
}
これらを組み合わせた最終形。
det :: [[Integer]] → Integer
det ass = head (head bss)
where
bss = foldl (trimult · mut) ass′ (replicate (n − 1) ass)
ass′ = if odd n then upper ass
else map (map negate) (upper ass)
n = length ass
Dart
// 方法3 の最終形： upper で上三角部分だけ保持し、trimult · mutTri を n-1 回反復
List<List<BigInt>> upper(List<List<BigInt>> ass) =>
[for (var i = 0; i < ass.length; i++) ass[i].sublist(i)];
BigInt detMatmultFinal(List<List<BigInt>> ass) {
final n = ass.length;
var bss = upper(ass);
if (n.isEven) {
bss = bss.map((r) => r.map((e) => -e).toList()).toList();
}
for (var i = 0; i < n - 1; i++) {
bss = trimult(mutTri(bss), ass);
}
return bss.first.first;
}
ここで upper = zipWith drop [0..] です。
Dart
// 章末: 3 通りの方法で同じ行列の行列式を計算してみる
void main() {
BigInt b(int x) => BigInt.from(x);
// 3x3 の例： |A| = -306（教科書的な例）
final a = [
[b(6), b(1), b(1)],
[b(4), b(-2), b(5)],
[b(2), b(8), b(7)],
];
print(detChio(a)); // -306
print(detChioImproved(a)); // -306
print(detMatmultFinal(a)); // -306
// 有理数版
Rat r(int x) => Rat(BigInt.from(x));
final ar = a.map((row) => row.map((e) => Rat(e)).toList()).toList();
final g = detGauss(ar);
print('\${g.n}/\${g.d}'); // -306/1
}
三つの方法の比較
要素が (−20, 20) の範囲にあるランダム行列を使って、いろいろな n で速度を比べた結果は次の通りです。
| 方法 | n = 150 の実行時間 |
|---|
| ガウスの消去法（有理数除算） | 約 30 秒 |
| 改良版キオ（凝縮＋除算を交互） | 約 10 秒 |
| 反復行列積（除算なし） | 約 40 秒 |
結論
予想どおり、元のキオ版はまったく使い物にならず（数字が巨大化する）、凝縮と除算を交互に行う改良版キオが明らかな勝者となりました。
結びに
アルゴリズムの出自
clow 列とは
置換のサイクル分解を一般化したもので、各サイクルが途中の要素の繰り返しを含んでもよい「閉じた歩行」のこと。Mahajan と Vinay は、置換に対応しない clow 列の符号付き項がすべて打ち消し合い、置換に対応する項 a1π(1)a2π(2) ··· anπ(n) だけが残ることを示した。この事実の純粋に代数的な証明はまだ知られていない。
配列を使っていないこと
本章のアルゴリズムは、Haskell の配列（可変・不変を問わず）をまったく使っていません。行列はいつも「行のリスト」でひっそりと表現されており、そのおかげで各アルゴリズムが簡潔に書けています。
ただし、より良い設計としては、行列に対する抽象データ型を定義して、「第 1 列」「第 1 行」「対角線」「主部分行列」などの操作をプリミティブとして提供する、というやり方も考えられます。
参考文献
Bareiss, E. H. (1968). Sylvester's identity and multi-step integer preserving Gaussian elimination. Mathematics of Computation 22 (103), 565–78.
Mahajan, M. and Vinay, V. (1997). Determinant: combinatorics, algorithms and complexity. Chicago Journal of Theoretical Computer Science, Article 5.