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

前回の記事ではMOLSを用いた直交表の構築方法について説明した。この方法では構築できる直交表のパターンに限界があるため、他のより優れた方法について調べてみた。

直交表の構築方法にはいろいろあるようだが、本稿では私が特に数学的に面白そうだと感じた有限射影幾何を用いた方法について説明したいと思う。ただし、有限射影幾何とはどういったものかという説明だけでもなかなかの文章量になりそうなので、やはり記事をいくつかに分けて説明したいと思う。その走りとして、本稿では手法の概要と有限射影幾何の基礎について述べた後、直交表の構築に必要な $G$ 行列というものを構築する方法について説明したいと思う。

$G$ 行列を用いた直交表の構築手法

手法の概要

論文[1]によると、強さ $t\ (t=2, 3, 4)\$ の直交表を作るためには、ある性質を持つ行列 $G$ があればよいらしい。その性質とは「任意の $t$ 行を選んだとき、それらが一次独立になっている」ということである。もしそのような行列 $G$ が得られれば、以下のようにして直交表が生成できるとのことだ。

$\displaystyle{ G\Theta }$

ただし、 $k$ を因子数とすると、 $G$ は $k$ 行 $d$ 列の行列、 $\Theta$ は $d$ 行 $q^d$ 列の行列である。特に、 $\Theta$ は $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間 $\mathbb{F}_q^d$ の相異なる元を列ベクトルとして $q^d$ 本並べることで得られる行列である。 $d$ は次元のパラメータであるが、詳細は後述する。

生成される行列 $G\Theta$ は $k$ 行 $q^d$ 列となる。この行列を転置すれば、前回の記事と同様のフォーマットの直交表が得られる。すなわち、 $(G\Theta)^t = \Theta^t G^t$ が求める直交表となる。

疑問点のオンパレードと解決の筋道

さて、これだけ見ると疑問点だらけでわけが分からない。少なくとも、以下のような疑問が湧いてくる。

$G$ 行列はどのようにして求めればよいのか？
$d$ というパラメータがいきなり出てきたが、これは何に基づいてどういう値に設定すべきものなのか？
なぜ $G \Theta$ が直交表になると言えるのか？
構築される直交表はどういうパターンになるか？

これらの疑問について、まずは現時点で私が理解していることの概要について述べる。

まず1点目について、これこそがまさに有限射影幾何が威力を発揮するところである。具体的には、有限体上の射影空間の点を求めて、それらから必要な部分集合を抜き出すことで解決できる。本稿の主題であり、難しいポイントなので詳細は後述する。

2点目について、dは因子数や水準数に依存して決まる。詳細は後述する。

3点目について、これは現時点で証明が見つけられていない。論文[1]で引用されている書籍に説明があるようだが、手元にないため参照できない。というわけで、今回も証明を自力で考えてみようと思うが、現時点での見通しはなく証明できるか分からない。少なくとも、これについては次回以降に回す。

4点目について、まず強さは2～4までいけるようである。5以上でも同様の理屈で行けそうだが、強さに応じて1点目で述べた「射影空間の点から必要な部分集合」なるものを抜き出すステップの計算量が爆発していくのと、実用上あまり必要になるとは思えないため、強さ5以上は考えないことにする。

有限体を利用している関係上、水準数は相変わらず素数冪に限られる。しかし、因子数は自由に選べるというのがこの手法の利点となる。

結局、構築できる直交表の型は $t\text{-}(q, k, \lambda)$ ( $2 \le t \le 4$ , $q$ は素数冪, $1 \le k, 1 \le \lambda)$ となる。 $\lambda$ は他のパラメータに応じて自動的に決まるため、無駄に大きくすることはできるが入力としての自由度は基本的にないものと考えてよい。

有限体上の射影空間

では、 $G$ 行列を求めるための有限射影幾何について考えていこう。具体的には有限体上の射影空間を考えるのだが、先に射影空間の一般論から説明する。

射影空間

射影空間の定義はいろいろなバリエーションがあるが、今回使用するものをWikipedia[2]から引用する。

Projective space

Given a vector space $V$ over a field $K$ , the projective space $P(V)$ is the set of equivalence classes of $V \setminus \{0\}$ under the equivalence relation $\sim$ defined by $x \sim y$ if there is a nonzero element $\lambda$ of $K$ such that $x = \lambda y$ .

つまり射影空間とは、互いにスカラー倍の関係にあるベクトル同士を同一視して、それらを改めて点と捉えた空間のことである。元のベクトル空間 $V$ の次元が $n+1$ なのに対して、そこから得られる射影空間の次元は $n$ であることに注意されたい。

また、射影空間における直線というものも考えることができる。まず、元の空間において原点を通る2つの異なる直線 $l_1, l_2$ があったとして、それらを両方とも含む平面 $S$ を考える。射影空間上では $l_1, l_2$ は点になるが、この平面 $S$ が射影空間上では点 $l_1, l_2$ を結ぶ直線となる。

もとのベクトル空間上での直線を点、平面を直線と捉えるため、少々ややこしい空間であるが、こういうのは慣れるしかない。

有限体の場合

さて、射影空間を得るためにはまずベクトル空間を考える必要があるが、体 $K$ としては何を選んでもよい。つまり、有限体上のベクトル空間であっても良いわけである。有限体上の射影空間のことを有限射影空間と呼ぶ。

以下では $\mathbb{F}_q$ ( $q$ は素数冪) 上の $n+1$ 次元ベクトル空間から得られる $n$ 次元射影空間を $\mathrm{PG}(n, q)$ と表記する。

有限射影空間の完全代表系

有限射影空間上の全ての点から完全代表系を求め、それを再び元のベクトル空間 $V$ の元と考えると、任意の2点から得られる $V$ 上のベクトルは一次独立となる。これは有限射影空間の定義を考えれば当然である。

有限射影空間の存在条件

有限射影空間は $q$ が素数冪の場合には存在することが知られているが、そうでない場合に存在するかどうかは未解決問題である[4]。そのため、以降の議論では $q$ が素数冪の場合のみを扱う。

有限射影空間上の点の数

$\mathbb{F}_q$ 上の $n$ 次元ベクトル空間 $V$ における点 (ベクトル) の数は $q^n$ 個に限られる。さらに、 $\mathrm{PG}(n, q)$ の場合は互いにスカラー倍の関係にある点は同一視されるので、さらに点の数が減ることになる。これは具体的にはいくつになるのだろうか？

次数が2の場合、つまり $\mathrm{PG}(2, q)$ 上の点の数は $q^2+q+1$ 個になるようである[3]。ちなみに、2次元の射影空間のことを特別に射影平面と呼ぶ。

次数が3以上の場合はどうなるのだろうか？これはGaussian binomial coefficientと呼ばれるものを使って以下のように算出できる[4][5]。

$\displaystyle{ \left( \begin{array}{c} n+1 \\ 1 \\ \end{array} \right)_q }$

私自身、Gaussian binomial coefficientなどというものは初めて見たのだが、計算の詳細は本稿の主眼ではないため、[5]などを参照されたい。ここではいくつかの具体的な $n, q$ について、 $\mathrm{PG}(n, q)$ 上の点の数がどうなるかをsageで計算してみよう。

sage: for dim in range(2,8):
....:     for order in range(2,8):
....:         if(is_prime_power(order)):
....:             print("The number of points in PG({}, {}) = {}".format(dim, order, q_binomial(dim+1, 1, q=order)))
....:
The number of points in PG(2, 2) = 7
The number of points in PG(2, 3) = 13
The number of points in PG(2, 4) = 21
The number of points in PG(2, 5) = 31
The number of points in PG(2, 7) = 57
The number of points in PG(3, 2) = 15
The number of points in PG(3, 3) = 40
The number of points in PG(3, 4) = 85
The number of points in PG(3, 5) = 156
The number of points in PG(3, 7) = 400
The number of points in PG(4, 2) = 31
The number of points in PG(4, 3) = 121
The number of points in PG(4, 4) = 341
The number of points in PG(4, 5) = 781
The number of points in PG(4, 7) = 2801
The number of points in PG(5, 2) = 63
The number of points in PG(5, 3) = 364
The number of points in PG(5, 4) = 1365
The number of points in PG(5, 5) = 3906
The number of points in PG(5, 7) = 19608
The number of points in PG(6, 2) = 127
The number of points in PG(6, 3) = 1093
The number of points in PG(6, 4) = 5461
The number of points in PG(6, 5) = 19531
The number of points in PG(6, 7) = 137257
The number of points in PG(7, 2) = 255
The number of points in PG(7, 3) = 3280
The number of points in PG(7, 4) = 21845
The number of points in PG(7, 5) = 97656
The number of points in PG(7, 7) = 960800

さらっと計算したが、sageは本当になんでも計算できてすごい。上記の結果から、 $n, k$ のいずれが増加しても点の数が爆発的に増加する様子が見て取れる。

有限射影空間と $G$ 行列

$G$ 行列の求め方

強さ2の場合

さて、話を直交表に戻そう。有限射影空間を用いることで、直ちに強さ2の直交表を求めることができる。というのも、前述の通り射影空間の完全代表系を元のベクトル空間に戻してやれば、そのうちの任意の2ベクトルは互いに一次独立となっており、これはまさに $G$ 行列の要件を満たすからである。

つまり、射影空間の完全代表系を元のベクトル空間に戻したものを行ベクトルとして縦に並べれば、それが $G$ 行列となる。

強さ3以上の場合

強さ3以上の場合は、まず強さ2の $G$ 行列を構築する。その上で、そこから任意の $t$ 行を選んだ時に、それらが一次独立になるように一生懸命探索して部分集合を選別するしかない。これをまじめにやると計算量が指数オーダーになってしまうので、以下のようなgreedyなアルゴリズムを使うなど、計算の工夫が必要となる。

選抜したベクトルの集合を空集合で初期化する。
強さ2の $G$ 行列の行を順に選んで以下を実行する。
今注目しているベクトルと現時点で選別されたベクトルの集合を合わせた集合について、そこから任意の $t$ 個のベクトルを選ぶと一次独立になるかを調べる。YESであれば今のベクトルを選別したベクトルの集合に加える。

ただし、これだと $G$ 行列の行数を最大化できないということが普通に起こりえる。また、それでなくとも完全代表系の選び方によっては $G$ 行列の行数を最大化できないケースがあるのではないのか？という疑問もある。しかし、実用上は $G$ の行数は多ければ多いほどよいが、必ずしも最大の行数となる必要はない。そのため、手法や計算の容易性のために $G$ 行列が多少小さくなることは許容できるものと考えている。

計算してみよう

射影空間上の点、及びその完全代表系はsageを使って計算が可能である。また、ベクトルの集合が一次独立であるかどうかもsageを使えば簡単に確かめられる。以下に $G$ 行列を求めるプログラムを示す。

#!/usr/bin/env sage

import sys
from sage.all import *

def getProjectiveSpacePoints(dim, order):
    PG = ProjectiveSpace(dim, GF(order, name="a", modulus="primitive"))
    return PG.rational_points()

def calcGMatrix(dim, order, strength):
    points = getProjectiveSpacePoints(dim, order)
    # Convert points in PG(dim, order) to vectors in vector space over F_{order}
    vectors = [vector(p) for p in points]
    if strength == 2:
        return Matrix(vectors)
    else:
        return Matrix(getLinearIndependentVectors(vectors, order, strength))

def getLinearIndependentVectors(vectors, order, strength):
    V = VectorSpace(GF(order, name="a", modulus="primitive"), len(vectors[0]))
    currentSet = []
    for v in vectors:
        if isRowIndependent(Matrix([v] + currentSet), order, strength):
            currentSet += [v]
    return currentSet

def genVectorIndex(index, numVecs):
    vectorsSet = []
    if len(index) == numVecs:
        vectorsSet.append([ind for ind in index])
    else:
        for i in range(len(index)-numVecs+1):
            if numVecs == 1:
                vectorsSet += [[index[i]]]
            else:
                subVectorsSet = genVectorIndex(index[i+1:], numVecs-1)
                for j in range(len(subVectorsSet)):
                    subVectorsSet[j] = [index[i]] + subVectorsSet[j]
                vectorsSet += subVectorsSet
    return vectorsSet

def isRowIndependent(G, order, strength):
    V = VectorSpace(GF(order, name="a", modulus="primitive"), len(G[0]))
    # If the number of rows of G matrix is relatively small,
    # all we have to do is just to check the linear independency of all rows of G.
    if G.nrows() <= strength:
        if V.linear_dependence([v for v in G]) != []:
            return False
    else:
        # If the dimension of the vector space is too small,
        # it is impossible to choose linear independent vector set.
        if G.ncols() < strength:
            return False
        vectorIndexSet = genVectorIndex([i for i in range(G.nrows())], strength)
        for vis in vectorIndexSet:
            vs = [G[i] for i in vis]
            if V.linear_dependence(vs) != []:
                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 = calcGMatrix(dim, order, strength)
    print("G matrix:")
    print(G)

    print("=== Sanity check ===")
    if isRowIndependent(G, order, strength):
        print("Correct G matrix.")
    else:
        print("Invalid G matrix.")

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

以下に実行例を示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py
usage: ./gmat.py dim order strength
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.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]
=== Sanity check ===
Correct G matrix.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 2 2 3
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
=== Sanity check ===
Correct G matrix.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 2 2 4
error: dim+1 should not be less than strength.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 3 3 4
G matrix:
[0 0 0 1]
[1 0 0 1]
[0 1 0 1]
[0 0 1 1]
[1 1 1 1]
=== Sanity check ===
Correct G matrix.

実行時引数として受け取っている次元の値は射影空間の次元としている。これは前述のパラメータ $d$ よりも1小さい値となっているので注意されたい。

得られた $G$ 行列は最後に任意の $t$ 行が一次独立になるかのチェックを行っているが、それは全てパスしている。

次元 $d$ の決め方

前述した疑問点の2つ目についてここで詳細に考えてみよう。直交表を作るにあたり、各因子は $G$ 行列の行に対応することになる。行を削る分には構わないが、行が足りないとなった場合には、次元を上げて行数を増やす必要がある。そのため、 $d, q$ の値をいろいろと変えたときに $G$ 行列の行数がどうなるかをあらかじめ調べておき、テスト要件として求められる因子数を満たすような値を指定すればよい。