直交表の数理 ～ MOLSを用いた直交表の構築編 ～

実験計画法

前回の記事ではmutually orthogonal Latin squares (MOLS) というテーブル群を構築する方法について説明した。本稿ではMOLSを用いて直交表を構築する方法について説明する。

直交表の構築

基本の構築方法

pを素数、 kを1以上の整数、 q = p^kとする。このとき、MOLS( q)からは 2\text{-}(q, m, 1) ( m = 2, \cdots, q+1) というタイプの直交表が得られる。以下ではMOLS(3)から直交表 2\text{-}(3, 4, 1)を得るケースを例にして、直交表の構築方法を説明する。全体的に英語版Wikipedia [1]を参考にした。

まず、直交表の最初の2列を順序対が辞書順になるように埋める。

# 因子1 因子2 因子3 因子4
1 1 1
2 1 2
3 1 3
4 2 1
5 2 2
6 2 3
7 3 1
8 3 2
9 3 3

次に、今埋めた2列の値をそれぞれLatin Squareの行・列と見なして、MOLSに属する各表から値を取り出す。今は例としてMOLS(3)を考えているので、これに属するLatin squareは L_1, L_2の2つである。そのため、3列目と4列目をそれぞれ L_1(i, j), L_2(i, j)の値で埋めればよい。

分かりやすさのために L_1, L_2を具体的に示しておく。

L_1

1 2 3
2 3 1
3 1 2

L_2

1 2 3
3 1 2
2 3 1

これより、以下のように直交表が得られる。

# 因子1 因子2 因子3 因子4
1 1 1 1 1
2 1 2 2 2
3 1 3 3 3
4 2 1 2 3
5 2 2 3 1
6 2 3 1 2
7 3 1 3 2
8 3 2 1 3
9 3 3 2 1

直交表が得られていることの証明

なんだかよくわからないうちに直交表が構築できてしまったように見えるのだが、これは本当に直交表になっているのだろうか？相変わらずWikipedia[1]には証明がないので、自分で考えてみた。なお、記述をシンプルにするため、これ以降はLatin Squareや直交表の行・列の番号は0始まりとする。

上記の方法によって得られるテーブルの行数は q^2となる。そのため、（強さ2で考えているので）任意の2列を選んで各行から順序対を形成したときに、それらがすべて異なることを示せばよい。そうすれば、ちょうど異なる q^2個の順序対が出尽くすことになるため、指数1の直交表が得られたことになる。

0列目、1列目、及び2列目以降で性質が異なるので、以下の組み合わせのパターンで場合分けをする。

  1. 0列目と1列目
  2. 0列目と2列目以降
  3. 1列目と2列目以降
  4. 2列目以降同士
準備

各ケースでの証明を始める前に、それぞれの列の i行目の値がどのように書き表せるかを最初にまとめておこう。有限体 \mathbb{F}_qの元は前回の記事に倣って a_0, \cdots, a_{q-1}と書く。

  • 0列目： \left\lfloor \frac{i}{q} \right\rfloor+1
  • 1列目： (i \bmod q)+1
  • r+1 ( 1 \le r \le q-1) 列目： a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor}
0列目と1列目

このケースは順序対が辞書順になるように値を決めているので自明である。

0列目と2列目以降

1列目と r+1 ( 1 \le r \le q-1) 列目について順序対を形成したとき、 i行目と i'行目の順序対が一致したと仮定する。すると、以下の2つの等式が成立する。

\displaystyle{ \begin{eqnarray} && \left\lfloor \frac{i}{q} \right\rfloor = \left\lfloor \frac{i'}{q} \right\rfloor \\ && a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} = a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }

1つ目の式を2つ目の式に代入して整理すると以下のようになる。

\displaystyle{ \begin{eqnarray} a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} \\ a_{(i \bmod q)} &=& a_{(i' \bmod q)} \\ \end{eqnarray} }

先ほどの1つ目の式より |i-i'| < qとなるため、 i \equiv i' \pmod qとなるためには i = i'となる他ない。

1列目と2列目以降

2列目と r+1 ( 1 \le r \le q-1) 列目について順序対を形成したとき、 i行目と i'行目の順序対が一致したと仮定する。すると、以下の2つの等式が成立する。

\displaystyle{ \begin{eqnarray} && (i \bmod q) = (i' \bmod q) \\ && a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} = a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }

1つ目の式を2つ目の式に代入して整理すると以下のようになる。

\displaystyle{ \begin{eqnarray} a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \\ a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }

2つ目の等号では r\ne 0より a_r \ne 0であることを用いた。

先ほどの1つ目の式より |i-i'|は qの倍数となるため、 \left\lfloor \frac{i}{q} \right\rfloor = \left\lfloor \frac{i'}{q} \right\rfloorとなるためには i = i'となる他ない。

2列目以降同士

2列目以降同士の組み合わせであれば、MOLSの定義より互いに直交するLatin squareから順序対を作ることになるため、それらは当然すべて異なる。

以上により、前述の方法で直交表が構築できることが分かった。

MOLSの限界とその他の手法

因子数の減らし方

前述した直交表の構築方法では、MOLSの完全集合の要素をすべて使用した。この場合、因子数が q+1に限定されてしまい、生成できる直交表のパターンがあまりにも限定的である。

実は、因子数を q+2以上に上げることは難しいが、下げることは割と簡単である。単にMOLSの要素から必要な分だけ排除してやれば良いだけである。例えば、先ほど直交表 2\text{-}(3, 4, 1)を構築したが、これを 2\text{-}(3, 3, 1)にしたければ単に L_2を無視して最終列を削れば良いだけである。

ただし、この場合でも行数を削ることはできない。これは、因子数が3だろうが4だろうが、例えば最初の2行の順序対のパターンを出し尽くすのに必要な行数は確保する必要があるためである。

MOLSの限界

結局、MOLS( q)から生成できる直交表のパターンは前述のとおり 2\text{-}(q, m, 1) ( m = 2, \cdots, q+1)となる。これがMOLSを用いた直交表構築法の限界である。具体的には、以下のような直交表は生成することができない。

  • 強さが3以上
  • 水準数が素数冪でない
  • 因子数が q+2以上

その他の直交表構築法

MOLSを用いる方法以外にも、直交表を構築する方法はあるようだ。詳しくは英語版Wikipedia[2]にいろいろ書いてある。例えば、Latin squareの代わりにLatin cubeとかLatin hypercubeと呼ばれるものを使う方法、アダマール行列なるものを使う方法、及び何やら難しい符号化の理論を使う方法があるらしい。いずれについても、現時点で私は理解できていない。

直交表を構築してみよう

では、実際に直交表を構築するプログラムを書いてみよう。前回の記事で作成したmols.pyを部品として使用し、直交表を構築するプログラム例を以下に示す。

#!/usr/bin/python3

import mols
import sys

def calcOrthogonalArray(ml, factor):
    order = len(ml[0])
    height = order**2
    oa = [[0 for i in range(factor)] for j in range(height)]
    for i in range(height):
        for j in range(factor):
            if j == 0:
                oa[i][j] = int(i/order) + 1
            elif j == 1:
                oa[i][j] = i%order + 1
            else:
                latSq = ml[j-2]
                oa[i][j] = latSq[int(i/order)][i%order]
    return oa

def main(argc, argv):
    if argc != 3:
        print("usage: " + argv[0] + " factor level")
        return

    factor = int(argv[1])
    level = int(argv[2])
    oa = []
    if factor == 1:
        # This is a special case
        oa = [[i+1] for i in range(level)]
    elif factor > level+1:
        print("error: factor should be less than or equal to level+1.")
    else:
        ml = mols.calcMOLS(level)
        oa = calcOrthogonalArray(ml, factor)

    for line in oa:
        print(line)

if __name__ == "__main__":
    main(len(sys.argv), sys.argv)

このプログラムの実行例を以下に示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 1 1
[1]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 1 3
[1]
[2]
[3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 2 3
[1, 1]
[1, 2]
[1, 3]
[2, 1]
[2, 2]
[2, 3]
[3, 1]
[3, 2]
[3, 3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 3 3
[1, 1, 1]
[1, 2, 2]
[1, 3, 3]
[2, 1, 2]
[2, 2, 3]
[2, 3, 1]
[3, 1, 3]
[3, 2, 1]
[3, 3, 2]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 4 3
[1, 1, 1, 1]
[1, 2, 2, 2]
[1, 3, 3, 3]
[2, 1, 2, 3]
[2, 2, 3, 1]
[2, 3, 1, 2]
[3, 1, 3, 2]
[3, 2, 1, 3]
[3, 3, 2, 1]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 5 3
error: factor should be less than or equal to level+1.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 6 5
[1, 1, 1, 1, 1, 1]
[1, 2, 2, 2, 2, 2]
[1, 3, 3, 3, 3, 3]
[1, 4, 4, 4, 4, 4]
[1, 5, 5, 5, 5, 5]
[2, 1, 2, 3, 4, 5]
[2, 2, 3, 5, 1, 4]
[2, 3, 5, 4, 2, 1]
[2, 4, 1, 2, 5, 3]
[2, 5, 4, 1, 3, 2]
[3, 1, 3, 4, 5, 2]
[3, 2, 5, 1, 4, 3]
[3, 3, 4, 2, 1, 5]
[3, 4, 2, 5, 3, 1]
[3, 5, 1, 3, 2, 4]
[4, 1, 4, 5, 2, 3]
[4, 2, 1, 4, 3, 5]
[4, 3, 2, 1, 5, 4]
[4, 4, 5, 3, 1, 2]
[4, 5, 3, 2, 4, 1]
[5, 1, 5, 2, 3, 4]
[5, 2, 4, 3, 5, 1]
[5, 3, 1, 5, 4, 2]
[5, 4, 3, 1, 2, 5]
[5, 5, 2, 4, 1, 3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 6 6
/usr/lib/python3/dist-packages/apport/report.py:13: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
  import fnmatch, glob, traceback, errno, sys, atexit, locale, imp, stat
Traceback (most recent call last):
  File "./oa.py", line 42, in <module>
    main(len(sys.argv), sys.argv)
  File "./oa.py", line 35, in main
    ml = mols.calcMOLS(level)
  File "/mnt/c/Users/shiny/Dropbox/programs/blog/2021/mols.py", line 17, in calcMOLS
    gf = GF(q, name="a", modulus="primitive")
  File "sage/structure/factory.pyx", line 367, in sage.structure.factory.UniqueFactory.__call__ (build/cythonized/sage/structure/factory.c:2162)
  File "/usr/lib/python3/dist-packages/sage/rings/finite_rings/finite_field_constructor.py", line 566, in create_key_and_extra_args
    raise ValueError("the order of a finite field must be a prime power")
ValueError: the order of a finite field must be a prime power

第1引数が因子数、第2引数が水準数である。因子数が q+1以下でなければならなかったり、水準数が素数冪でなければならない制限が見て取れる*1。

まとめ

本稿ではMOLSから直交表を構築する方法について説明した。また、この方法で生成できる直交表のバリエーションの限界について述べた。

現時点では良く分かっていないのだが、MOLSは有限射影幾何なるものと関連があるらしい。次回はそれについて調べてまとめてみたいと思う。ただ、見通しはないのでどうなるかは分からない。調べた結果、私が理解できて、かつ面白いと思ったら記事を書こうと思う。

参考

[1] Mutually orthogonal Latin squares - Wikipedia
[2] Orthogonal array - Wikipedia

*1:水準数が素数冪でない場合のエラーハンドリングがなく、sageの関数内で異常終了しているが、サンプルプログラムなのでご容赦願いたい。