2^13次正方行列の積の計算は数十秒でできる

例のSNS上で、こんな投稿を見つけた。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">n=2^13次正方行列の積ってGPUしても終わらないよね？現実的な時間で終わった人いますか？ <a href="https://t.co/AglBRaDQcY">pic.twitter.com/AglBRaDQcY</a></p>&mdash; りーさぶさぶ (@li__subsub) <a href="https://twitter.com/li__subsub/status/1757996023983771769?ref_src=twsrc%5Etfw">February 15, 2024</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

「n=2^13次正方行列の積」が「現実的な時間」で終わったかどうかを聞いているようである。

# 定義の確認

## n=2^13次正方行列

「n=2^13次正方行列」とは何かを確認する。

> 正方行列（せいほうぎょうれつ、英: square matrix）とは、行要素の数と列要素の数が一致する行列である。サイズが n × n つまり、n 行 n 列であるとき、n 次正方行列という。 

[正方行列 - Wikipedia](https://ja.wikipedia.org/wiki/%E6%AD%A3%E6%96%B9%E8%A1%8C%E5%88%97)

「2^13」に用いられている演算子 `^` は、C言語などでは xor (排他的論理和) を表す。
しかし、最初に紹介した投稿には Wolfram Alpha のスクリーンショットが貼られており、「2^13」は 2 の 13 乗 $2^{13}$ を表すと考えられる。
$2^{13} = 8192$ である。

よって、「n=2^13次正方行列」とは、8192 行 8192 列の行列であるといえる。

## 積

> 積（せき）とは数学の乗法の結果を指す。

[積 - Wikipedia](https://ja.wikipedia.org/wiki/%E7%A9%8D)

「積」とは、掛け算の答えである。
最初に紹介した投稿に行列の行数 (=列数) の 3 乗の計算を行っている様子が載っていることから、今回は[行列の積](https://ja.wikipedia.org/wiki/%E8%A1%8C%E5%88%97%E3%81%AE%E4%B9%97%E6%B3%95)の計算を表すと考えられる。
この行列の積は、$i$ 行 $k$ 列の行列と $k$ 行 $j$ 列の行列から $i$ 行 $j$ 列の行列を得る計算であり、定義どおりに素直に計算した際の計算量は $O(ijk)$ になる。

## 現実的な時間

「現実的な時間」の定義はよくわからない。
適当にググったところ、以下の記述を見つけた。

> 8桁の英小文字+数字の組み合わせでも2分13秒という現実的な時間で解読できた。 

[6桁の英字パスワードなら家庭用PCでも1秒で解読。PPAPに警鐘 - PC Watch](https://pc.watch.impress.co.jp/docs/news/1333357.html)

よって、少なくとも2分13秒以内であれば「現実的な時間」であるといえるだろう。

# 計算量の見積もり

最初に紹介した投稿にも載っている通り、$n=2^{13}$ のとき $O(n^3)$ の計算を行うと、そのコストはだいたい $5.5 \times 10^{11}$ となる。

> AtCoderの問題では実行時間制約が2秒くらいであることが多いです。 コンテストに参加する人は、1秒あたり$10^8$回くらいの計算ができることを覚えておきましょう。

[W - 2.06.計算量](https://atcoder.jp/contests/apg4b/tasks/APG4b_w)

有名な競技プログラミングサイトの AtCoder では、1秒あたり $10^8$ 程度のコストの計算ができるというのが目安になっている。
ここでは、プログラムを**シングルスレッド**で実行し、**数秒以内**に実行を完了させることを目指す。
$5.5 \times 10^{11}$ を $10^8$ で割ると約 $5497$ となり、これは計算完了までに**約1.5時間**かかることが予測されることを示している。
これでは、通常の競技プログラミングでは実用的ではないだろう。

では、競技プログラミングではなかったらどうだろうか。
CTF（Capture The Flag）では、問題のデータを手元で解析し、求まった答え (flag) だけをスコアサーバに送ればいい問題がある。
この場合、数秒で実行が完了する必要は無く、競技時間内に計算を完了させて正しい答えを提出できれば得点が得られる。
競技時間は24時間や48時間に設定されていることが多いため、数時間程度で計算が完了するのであれば**十分実用的である**といえる。
さらに、このときの解析には**マルチスレッド**も利用することができる。
たとえば、60並列で計算を行うと、シングルスレッドで1.5時間かかる計算が**1.5分**で終わることが期待できる。
(各種オーバーヘッドや並列化しにくい部分の影響により、このような理想的な時間短縮ができない可能性もある)

よって、CTF の文脈では、$2^{13} \times 2^{13}$ の行列同士の積の計算は**十分実用的な時間で完了させることができる**と予測できる。

# 実験

実際に行列の積の計算を行うプログラムを実装し、実験を行ってみた。
「行列」の中身や計算の回数の情報は最初に紹介した投稿には無いようなので、今回は以下の条件で計算を行う。

* 行列の要素は `float` 型
* 計算前の行列の各要素の値は 0～1 の範囲
* 積の計算は1回だけ行う

以下が、今回実装したプログラムである。

```c:matrix_mult_test.c
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char* argv[]) {
	int n;
	int i;
	float *a, *b, *prod;
	size_t matrix_size, si;
	int adler_a, adler_b;
	if (argc < 2 || sscanf(argv[1], "%d", &n) != 1 || n <= 0) {
		fprintf(stderr, "Usage: %s n\n", argc > 0 ? argv[0] : "matrix_mult_test");
		return 1;
	}
	/* メモリを確保する */
	matrix_size = sizeof(float) * n * n;
	a = malloc(matrix_size);
	b = malloc(matrix_size);
	prod = malloc(matrix_size);
	if (a == NULL || b == NULL || prod == NULL) {
		fputs("memory allocation failed\n", stderr);
		free(a);
		free(b);
		free(prod);
		return 2;
	}
	/* 行列を適当に初期化する */
	for (i = 0; i < n; i++) {
		int j;
		for (j = 0; j < n; j++) {
			a[i * n + j] = (float)rand() / RAND_MAX;
			b[i * n + j] = (float)rand() / RAND_MAX;
		}
	}
	/* 行列の積の計算を実行する */
	for (i = 0; i < n; i++) {
		int j;
		for (j = 0; j < n; j++) {
			prod[i * n + j] = 0;
		}
	}
	#ifdef _OPENMP
	#pragma omp parallel for
	#endif
	for (i = 0; i < n; i++) {
		int j, k;
		for (k = 0; k < n; k++) {
			for (j = 0; j < n; j++) {
				prod[i * n + j] += a[i * n + k] * b[k * n + j];
			}
		}
	}
	/* 計算結果の Adler-32 を求める (最適化避け) */
	adler_a = 1;
	adler_b = 0;
	for (si = 0; si < matrix_size; si++) {
		adler_a = (adler_a + ((unsigned char*)prod)[si]) % 65521;
		adler_b = (adler_b + adler_a) % 65521;
	}
	printf("%04x%04x\n", adler_b, adler_a);
	free(a);
	free(b);
	free(prod);
	return 0;
}
```

行列の積の計算の高速化には、以下のテクニックなどが知られている。

* 小さなブロックごとに計算を行う
* SIMD命令を活用する
* データのプリフェッチを行う

しかし、今回はこのような複雑なテクニックの実装は避け、以下の単純なテクニックのみを取り入れた。

* OpenMP による並列化を行う
* キャッシュを考慮し、メモリ上で隣に配置された要素を連続でアクセスしていくループの順序にする

以下の手元の環境で実行した。

* CPU: 13th Gen Intel(R) Core(TM) i7-13700H 2.40GHz
  * コア数: 14、論理プロセッサ数: 20
* RAM: 32.0 GB
* OS: Windows 11 Home 22H2

コンパイラのバージョンとオプションは以下である。

```
YUKI.N>clang --version
clang version 17.0.6 (https://github.com/llvm/llvm-project.git 6009708b4367171ccdbf4b5905cb6a803753fe18)
Target: x86_64-w64-windows-gnu
Thread model: posix
InstalledDir: C:/MyApps/llvm-mingw-20231128-msvcrt-x86_64/bin

YUKI.N>clang -O3 -march=native -o matrix_mult_test.exe matrix_mult_test.c -Wall -Wextra -fopenmp
```

環境変数による OpenMP の動作設定は行っていない。

計算のサイズ (第1引数) を `8192` に設定し、[ストップウォッチD](https://www.vector.co.jp/soft/winnt/personal/se409211.html) を用いて実行にかかった時間を測定した結果、**約33秒** で実行が完了した。
この時間には行列の値の設定やチェックサムの計算などの積の計算以外の処理にかかる時間も含まれるため、積の計算は**30秒程度**でできたといえるだろう。

# 結論

最近のパソコンには、$2^{13}$ (8192) 次正方行列同士の積の計算を1分以内にできるものがある。