直交表の数理～有限射影幾何を用いた直交表の構築編～

前回の記事で有限射影空間の理論を用いた $G$ 行列の構築について説明した。本稿では前回説明した手法で正しく直交表が構築できることを証明し、実際に直交表を構築するプログラム例を示したいと思う。

有限射影幾何を用いた直交表構築手法の証明

おさらい

$k$ を因子数、 $t$ を直交表の強さ、 $q$ を素数冪とする。 $k$ 行 $d$ 列の行列 $G$ があって、その任意の $t$ ( $t \le d$ ) 行は $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間上で一次独立であるとする。このとき、 $d$ 行 $q^d$ 列の行列 $\Theta$ を用いて、直交表は $G\Theta$ という行列積により得られるとのことだった（詳細は前回の記事を参照）。以下ではこの事実の証明を考えていく。

事前に断らせて頂くと、証明はほぼ全て私が考えたものだが、有限体上のベクトル空間の扱いにあまり自信がない。有限体上のベクトル空間では、例えば内積が0でもベクトルが直交しないようなケースがあるなど、標数0の体で使える議論が使えないケースがある[1]。もしそういった点で誤りがあれば、Twitterやブログコメントにてご指摘頂けるとありがたい。

証明

${\bf x}_1, \cdots , {\bf x}_k$ を $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間 $V$ 上の行ベクトルの集合とし、これらのうち任意の $t$ 個を選べばそれらは一次独立になっているとする。これを用いて $G$ は以下のように書ける。

$\displaystyle{ G = \left( \begin{array}{cc} {\bf x}_1 \\ \vdots \\ {\bf x}_k \end{array} \right) }$

同様に ${\bf y}_1, \cdots , {\bf y}_{q^d}$ を $V$ 上の列ベクトルの集合とし、これらは相異なるとする。ベクトルの数が $q^d$ であるため、これらのベクトルは $V$ のベクトル全体の集合となる。これを用いて $\Theta$ は以下のように書ける。

$\displaystyle{ \Theta = \left( \begin{array}{ccc} {\bf y}_1 & \cdots & {\bf y}_{q^d} \end{array} \right) }$

これより、行列積 $G \Theta$ は以下のようになる。

$\displaystyle{ G \Theta = \left( \begin{array}{ccc} {\bf x}_1 {\bf y}_1 & \cdots & {\bf x}_1 {\bf y}_{q^d} \\ \vdots & \ddots & \vdots \\ {\bf x}_k {\bf y}_1 & \cdots & {\bf x}_k {\bf y}_{q^d} \end{array} \right) }$

示したいことは、 $G \Theta$ の任意の $t$ 行を選んだとき、それらの各列の値から順序集合を形成したとき、同じ順序集合がちょうど $q^{d-t}$ 回ずつ現れることである。というのも、列は全部で $q^d$ 個あり、可能な順序集合のパターンは全部で $q^t$ 通りとなるため、 $t$ 個の値の組がそれぞれ $q^d/q^t = q^{d-t}$ 回ずつ出現すれば直交表としての性質を満たすためである。

以下、選ばれた $t$ 個の行の番号を $i_1, \cdots , i_t$ とする。

まず、全ての順序集合が $q^{d-t}$ 回未満しか出現しないと仮定する。このとき、可能な順序集合の組み合わせが $q^t$ 通りしかないため、最大でも $q^t (q^{d-t}-1) = q^d - q^t$ 個の列しか存在できなくなってしまう。そのため仮定は誤りで、少なくとも1つの順序集合は $q^{d-t}$ 回以上出現する。

一番出現回数が多い順序集合について、それが出現する列がちょうど $q^{d-t}$ 個であることを示す。そのために、このような列が $q^{d-t}+1$ 個あると仮定して矛盾を導く。今、 $l_1, \cdots , l_{q^{d-t}+1}$ という番号の列について、 $i_1, \cdots , i_t$ 行目の値から成る順序集合が一致したと仮定する。この時、 ${\bf x}_{i_1}$ について以下の式が成り立つ。

$\displaystyle{ \begin{eqnarray} {\bf x}_{i_1} {\bf y}_{l_\alpha} &=& {\bf x}_{i_1} {\bf y}_{l_\beta} \\ {\bf x}_{i_1} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \end{eqnarray} }$

ただし、 $\alpha , \beta = 1, \cdots , q^{d-t}+1$ である。

${\bf x}_{i_2}, \cdots , {\bf x}_{i_t}$ についても同様に以下の式が得られる。

$\displaystyle{ \begin{eqnarray} {\bf x}_{i_2} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \\ &\vdots& \\ {\bf x}_{i_t} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \end{eqnarray} }$

この関係を線形写像の観点から整理してみる。そのために以下のような行列 $A$ を考える。

$\displaystyle{ A = \left( \begin{array}{c} {\bf x}_{i_1} \\ \vdots \\ {\bf x}_{i_t} \end{array} \right) }$

${\bf x}_{i_1}, \cdots , {\bf x}_{i_t}$ は一次独立なので、 $\mathrm{rank}A = t$ である。そのため、次元定理より $\mathrm{ker}A = d-t$ となる。 ${\bf y}_{l_\alpha} - {\bf y}_{l_\beta} \in \mathrm{ker}A$ なので、 ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ は $\mathbb{F}_q^d$ において全て同一の $d-t$ 次元空間に含まれる。そこで、 $\mathrm{ker}A$ のベクトル空間としての基底を ${\bf z}_1, \cdots , {\bf z}_{d-t}$ とすると、 ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ は全て以下の形式で表すことができる。

${\bf y}_{l_1} = c_1 {\bf z}_1 + \cdots + c_{d-t} {\bf z}_{d-t}$

ただし、 $c_1, \cdots , c_{d-t} \in \mathbb{F}_q$ である。これらの係数は $0, \cdots, q-1$ の $q$ 通りの値を取るため、この形式で表されるベクトルは全部で $q^{d-t}$ 本となる。これは ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ が相異なる $q^{d-t}+1$ 個のベクトルであることに矛盾する。

よって仮定は誤りで、一番出現回数が多い列は $q^{d-t}$ 個または $q^{d-t}+2$ 個以上ということになる。後者は $q^{d-t}+1$ 個の場合と同様の理由でありえないため、結局列数は $q^{d-t}$ 個ということになる。

今考えた一番出現回数が多い順序集合が出現する列を除いて、他の列について同様の議論を繰り返すことで、結局全ての順序集合の出現回数は $q^{d-t}$ 回と言える。よって $G \Theta$ は直交表になっている。

計算してみよう

証明ができて気分がすっきりしたところで、実際に計算してみよう。前回の記事で作成したgmat.pyをimportして直交表を構築するプログラム例を以下に示す。

#!/usr/bin/env sage

import gmat
import sys
from sage.all import *

def calcThetaMatrix(dim, order):
    K = GF(order, name="a", modulus="primitive")
    gamma = K.gen()
    Theta = calcThetaMatrixHelper(dim, order, gamma, 0, [0 for d in range(dim)])
    Theta = Matrix(Theta)
    return Theta.transpose()

def calcThetaMatrixHelper(dim, order, gamma, pos, currentVec):
    if pos == dim:
        return [currentVec.copy()]

    currentVec[pos] = 0
    Theta = calcThetaMatrixHelper(dim, order, gamma, pos+1, currentVec)
    for q in range(1, order):
        currentVec[pos] = gamma**q
        Theta += calcThetaMatrixHelper(dim, order, gamma, pos+1, currentVec)

    return Theta

def convertToInteger(oa, order):
    gf = GF(order, name="a", modulus="primitive")
    gamma = gf.gen()
    orderedGf = [0] + [gamma**(val) for val in range(order-1)]
    # Hasmap to convert elements of GF(q) to corresponding integers
    gfToInt = {}
    for i, val in enumerate(orderedGf):
        gfToInt[val] = i+1

    intOa = []
    for i in range(oa.nrows()):
        row = []
        for j in range(oa.ncols()):
            row.append(gfToInt[oa[i][j]])
        intOa.append(row)
    return Matrix(intOa)

def isCorrectOA(oa, strength):
    vectorIndexSet = gmat.genVectorIndex([i for i in range(len(oa[0]))], strength)
    for vis in vectorIndexSet:
        tupleCount = {}
        for row in range(oa.nrows()):
            rowTuple = tuple([oa[row][val] for val in vis])
            if rowTuple in tupleCount:
                tupleCount[rowTuple] += 1
            else:
                tupleCount[rowTuple] = 1

        count = 0
        for v in tupleCount.values():
            if count == 0:
                count = v
            elif count != v:
                return False
    return True

def main(argc, argv):
    if argc != 4:
        print("usage: " + argv[0] + " dim order strength")
        return

    dim = int(argv[1])
    order = int(argv[2])
    strength = int(argv[3])
    if order < 2 or dim < 1 or strength < 2:
        print("error: invalid parameter:")
        print("    dim={}, order={}, strength={}".format(dim, order, strength))
        return
    if dim + 1 < strength:
        print("error: dim+1 should not be less than strength.")
        return

    G = gmat.calcGMatrix(dim, order, strength)
    print("G matrix:")
    print(G)
    Theta = calcThetaMatrix(dim+1, order)
    print("Theta matrix:")
    print(Theta)
    oa = G * Theta
    oa = oa.transpose()
    intOa = convertToInteger(oa, order)
    intOa = Matrix(sorted(intOa))
    print("Resulting OA:")
    print(intOa)

    print("=== Sanity check ===")
    if isCorrectOA(intOa, strength):
        print("Correct OA.")
    else:
        print("Invalid OA.")

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

$\Theta$ を求めるところ（calcThetaMatrix関数）がややこしいが、単純なbacktrackアルゴリズムで生成している。

以下に実行例を示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py
usage: ./pg.py dim order strength
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 2 2
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
[0 1 0]
[1 1 0]
[1 0 0]
Theta matrix:
[0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1 1 1 1]
[1 1 2 2 2 2 1]
[1 2 1 2 1 2 2]
[1 2 2 1 2 1 2]
[2 1 1 2 2 1 2]
[2 1 2 1 1 2 2]
[2 2 1 1 2 2 1]
[2 2 2 2 1 1 1]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 2 3
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
Theta matrix:
[0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1]
[1 1 2 2]
[1 2 1 2]
[1 2 2 1]
[2 1 1 2]
[2 1 2 1]
[2 2 1 1]
[2 2 2 2]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 3 2 2
G matrix:
[0 0 0 1]
[1 0 0 1]
[0 1 0 1]
[1 1 0 1]
[0 0 1 1]
[1 0 1 1]
[0 1 1 1]
[1 1 1 1]
[0 0 1 0]
[1 0 1 0]
[0 1 1 0]
[1 1 1 0]
[0 1 0 0]
[1 1 0 0]
[1 0 0 0]
Theta matrix:
[0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1]
[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[1 1 1 1 2 2 2 2 2 2 2 2 1 1 1]
[1 1 2 2 1 1 2 2 1 1 2 2 2 2 1]
[1 1 2 2 2 2 1 1 2 2 1 1 2 2 1]
[1 2 1 2 1 2 1 2 1 2 1 2 1 2 2]
[1 2 1 2 2 1 2 1 2 1 2 1 1 2 2]
[1 2 2 1 1 2 2 1 1 2 2 1 2 1 2]
[1 2 2 1 2 1 1 2 2 1 1 2 2 1 2]
[2 1 1 2 1 2 2 1 2 1 1 2 2 1 2]
[2 1 1 2 2 1 1 2 1 2 2 1 2 1 2]
[2 1 2 1 1 2 1 2 2 1 2 1 1 2 2]
[2 1 2 1 2 1 2 1 1 2 1 2 1 2 2]
[2 2 1 1 1 1 2 2 2 2 1 1 2 2 1]
[2 2 1 1 2 2 1 1 1 1 2 2 2 2 1]
[2 2 2 2 1 1 1 1 2 2 2 2 1 1 1]
[2 2 2 2 2 2 2 2 1 1 1 1 1 1 1]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 4 3
G matrix:
[    0     0     1]
[    a     0     1]
[    0     a     1]
[    a     a     1]
[    a     1     0]
[a + 1     1     0]
Theta matrix:
[    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     a     a     a     a     a     a     a     a     a     a     a     a     a     a     a     a a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1]
[    0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1]
[    0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1]
Resulting OA:
[1 1 1 1 1 1]
[1 1 2 2 4 4]
[1 1 3 3 2 2]
[1 1 4 4 3 3]
[1 2 1 2 2 3]
[1 2 2 1 3 2]
[1 2 3 4 1 4]
[1 2 4 3 4 1]
[1 3 1 3 3 4]
[1 3 2 4 2 1]
[1 3 3 1 4 3]
[1 3 4 2 1 2]
[1 4 1 4 4 2]
[1 4 2 3 1 3]
[1 4 3 2 3 1]
[1 4 4 1 2 4]
[2 1 1 2 3 2]
[2 1 2 1 2 3]
[2 1 3 4 4 1]
[2 1 4 3 1 4]
[2 2 1 1 4 4]
[2 2 2 2 1 1]
[2 2 3 3 3 3]
[2 2 4 4 2 2]
[2 3 1 4 1 3]
[2 3 2 3 4 2]
[2 3 3 2 2 4]
[2 3 4 1 3 1]
[2 4 1 3 2 1]
[2 4 2 4 3 4]
[2 4 3 1 1 2]
[2 4 4 2 4 3]
[3 1 1 3 4 3]
[3 1 2 4 1 2]
[3 1 3 1 3 4]
[3 1 4 2 2 1]
[3 2 1 4 3 1]
[3 2 2 3 2 4]
[3 2 3 2 4 2]
[3 2 4 1 1 3]
[3 3 1 1 2 2]
[3 3 2 2 3 3]
[3 3 3 3 1 1]
[3 3 4 4 4 4]
[3 4 1 2 1 4]
[3 4 2 1 4 1]
[3 4 3 4 2 3]
[3 4 4 3 3 2]
[4 1 1 4 2 4]
[4 1 2 3 3 1]
[4 1 3 2 1 3]
[4 1 4 1 4 2]
[4 2 1 3 1 2]
[4 2 2 4 4 3]
[4 2 3 1 2 1]
[4 2 4 2 3 4]
[4 3 1 2 4 1]
[4 3 2 1 1 4]
[4 3 3 4 3 2]
[4 3 4 3 2 3]
[4 4 1 1 3 3]
[4 4 2 2 2 2]
[4 4 3 3 4 4]
[4 4 4 4 1 1]
=== Sanity check ===
Correct OA.

上記のプログラムでは $G \Theta$ を転置して $\Theta^t G^t$ の計算結果を出力している点、及び入力する次元の値は射影空間の次元の値としている (つまり $d-1$ の値を入力する仕様である) 点に注意されたい。

また、有限体の元をそのまま出力すると分かりづらいので、以前の記事と同じ方法で整数値に変換している。その後、行を昇順にsortしたものを最終出力としている。

今回も最後に正しく直交表が生成できているかチェックしているが、全てパスしている。

ちなみに

最近気づいたのだが、実はsageには直交表を直接的に求めるための関数が存在している[2]。「なんだ、こんな便利な関数があるのか」とガックリしたのだが、少し触ってみたところ、指数2以上の直交表は構築できないような挙動をしている。単なる私の理解不足かもしれないが、なんだかいまいちな仕様であり、本稿で作ったプログラムも無意味ではなかったかなと思っている。

以下に実行例を示す。第1引数が因子数、第2引数が水準数で、 $t$ というパラメータに強さの値を与えればよい。

sage: for row in designs.orthogonal_arrays.build(2, 2, t=2):
....:     print(row)
....:
[0, 0]
[0, 1]
[1, 0]
[1, 1]
sage: for row in designs.orthogonal_arrays.build(3, 2, t=2):
....:     print(row)
....:
[0, 0, 0]
[0, 1, 1]
[1, 0, 1]
[1, 1, 0]
sage: for row in designs.orthogonal_arrays.build(4, 2, t=2):
....:     print(row)
....:
---------------------------------------------------------------------------
EmptySetError                             Traceback (most recent call last)
<ipython-input-156-349912d12039> in <module>
----> 1 for row in designs.orthogonal_arrays.build(Integer(4), Integer(2), t=Integer(2)):
      2     print(row)
      3

/usr/lib/python3/dist-packages/sage/combinat/designs/orthogonal_arrays.py in build(k, n, t, resolvable)
   2142
   2143         """
-> 2144         return orthogonal_array(k,n,t,resolvable=resolvable)
   2145
   2146     @staticmethod

/usr/lib/python3/dist-packages/sage/combinat/designs/orthogonal_arrays.py in orthogonal_array(k, n, t, resolvable, check, existence, explain_construction)
    891         if explain_construction:
    892             return msg
--> 893         raise EmptySetError(msg)
    894
    895     elif k <= t:

EmptySetError: There exists no OA(4,2) as k(=4)>n+t-1=3