三角行列のインデックス - ロベールの小部屋

　巨大な行列を扱っていると、メモリ容量を軽減するために色々と工夫を行う必要が出てくる。自分の専門であるところの量子化学計算もまさにそういうものであり、特に一番活用されるのが三角行列の扱いである。
　三角行列とは正方行列の対角線を境界として、どちらか一方が全て０であるような行列のことである。対角と右上のみに要素がある場合は上三角行列、対角と左下のみに要素がある場合は下三角行列と言う。実際には純粋な三角行列を扱うだけではなく、対称行列や反対称行列などの重複部分を無視し、残りの部分を三角行列と見做して利用する場合もある。
　三角行列を一次元配列にパッキングする際のインデクシングには大まかに２つ存在する。１つは、上三角行列の左上から右向き優先で順番に番号を振っていくタイプ。もう１つは下三角行列の左上から右向き優先で順番に番号を振っていくタイプである。例を言うならば、前者の場合、 $A_{00}$ , $A_{01}$ , $A_{02}$ , ..., $A_{0N}$ , $A_{11}$ , $A_{12}$ , ..., $A_{1N}$ , ..., $A_{NN}$ という順番で通し番号を振っていくものであり、後者の場合は、 $A_{00}$ , $A_{10}$ , $A_{11}$ , $A_{20}$ , $A_{21}$ , $A_{22}$ , ..., $A_{N0}$ , $A_{N1}$ , ..., $A_{NN}$ という順番で通し番号を振っていくものである（ $A_{ij}$ において、 $i$ は行番号、 $j$ は列番号とする）。それぞれ下向き優先に番号を振っていく事もできるが、それは単に転置（ここでは上三角行列と下三角行列の相互変換に相当）してから右向き優先で番号を振ったものと考えれば良いので、特に新しいインデクシングというわけではない。
　では、それぞれのインデクシングで、具体的に通し番号はどうなるのだろうか？　このインデックスを得るには「それまでに何個の要素があったか」を数えれば良い。
　先ず、具体的な計算の前に、いくつかの定義を行う。

$N+1$ 次元上三角行列 $A$ の要素 $A_{ij}$ のインデックスを $U(i, j: N)$ と表記する。
$N+1$ 次元下三角行列 $A$ の要素 $A_{ij}$ のインデックスを $D(i, j: N)$ と表記する。
$i$ , $j$ , $U(i, j: N)$ , $D(i, j: N)$ は全て 0 オリジンとする。

　では、上三角行列において計算してみよう。先ず、 $i$ 番目の行の自明に非 0 となる要素は $A_{ii}$ から $A_{iN}$ までの $N-i+1$ 個である。従って、 $A_{ii}$ までの要素数は 0 行目から $i-1$ 行目までの要素数を総計すれば簡単に求められる。

　 $U(i, i: N)$
　　 $= \sum_{n=0}^{i-1} (N-n+1)$
　　 $= (N+1)i - \frac{i(i-1)}{2}$
　　 $= \frac{i(2N+3-i)}{2}$

ここまで求まれば、 $U(i, j: N)$ は $U(i, i: N)$ から $j-i$ 進んだだけなので、

　 $U(i, j: N)$
　　 $= U(i, i: N) + j-i$
　　 $= \frac{i(2N+1-i)}{2} + j$

と求める事ができる。
　次に、下三角行列でのインデックスを求める。こちらの場合、 $i$ 番目の行の要素数は $i+1$ 個なので、 $A_{i0}$ までの要素数は

　 $D(i, 0: N)$
　　 $= \sum_{n=0}^{i-1} (i+1)$
　　 $= \frac{i(i-1)}{2} + i$
　　 $= \frac{i(i+1)}{2}$

となる。 $D(i, j: N)$ は $D(i, 0: N)$ から $j$ 進んだだけなので、ただちに

　 $D(i, j: N)$
　　 $= D(i, 0: N) + j$
　　 $= \frac{i(i+1)}{2} + j$

と求まる。
　 $U(i, j: N)$ と $D(i, j: N)$ を見比べてみると、ある事に気づくだろう。 $D(i, j: N)$ は、表記に $N$ という変数が書いてはいるものの、全く $N$ に依存していないのである。これは下三角行列のインデクシングの順番を見ても一目瞭然であり、行列のサイズに無関係に利用する事ができるようになっている。また、式も下三角行列の方が簡単になっている。この特徴から、実用的には専ら下三角行列タイプのインデクシングが利用される。ただ、状況によっては上三角行列の方が扱いやすくなる可能性はある（上三角行列の形式だと既存のあるファイルからリニアに読み込める状況など）。

　では、添字が３つ、あるいはもっと一般的に添字が $M$ 個の行列ではどうなるだろうか？　この事を考える前に、もう少し定義を行う。

添字が $M$ 個の $N+1$ 次元超三角行列 $A$ の要素 $A_{ij...k}$ のインデックスを $U(i, j, ..., k: N, M)$ あるいは $D(i, j, ..., k: M)$ とする。超下三角行列の表式は $N$ に依存しないため、 $N$ を書かないこととする。
超上三角行列は $A_{00...00}$ , $A_{00...01}$ , $A_{00...02}$ , ..., $A_{00...0N}$ , $A_{00...11}$ , ..., $A_{00...1N}$ , ..., $A_{00...NN}$ , ..., $A_{0N...NN}$ , $A_{11...11}$ , ..., $A_{11...1N}$ , ..., $A_{1N...NN}$ , ..., $A_{22...22}$ , ..., $A_{NN...NN}$ の順番で通し番号がふられる。

　つまり、添字は常に右の方が大きい（または等しい）状態で、右の方の添字から優先的に値を大きくしていく。

超下三角行列は $A_{000...0}$ , $A_{100...0}$ , $A_{110...0}$ , $A_{111...0}$ , ..., $A_{111...1}$ , $A_{200...0}$ , $A_{210...0}$ , $A_{211...0}$ , ..., $A_{211...1}$ , $A_{220...0}$ , ..., $A_{222...2}$ , ..., $A_{NN...NN}$ の順番で通し番号がふられる。

　つまり、添字は常に左の方が大きい（または等しい）状態で、一番右まで左から順番に値を増やしていき、端に到達したら左の方の値を増やして、以降を 0 にリセットして、また同様に増やしていく（分かるかな？）。
　超下三角行列の特徴は、0 から $n$ のみの添字を使ったパターンは $A_{nnn...n}$ までに網羅されるということである。

　では、超上三角行列から考えてみよう。例えば $N = 2$ , $M = 3$ の場合、

　 $A_{000}$ , $A_{001}$ , $A_{002}$ , $A_{011}$ , $A_{012}$ , $A_{022}$ , $A_{111}$ , $A_{112}$ , $A_{122}$ , $A_{222}$

の順に番号が振られることになる。ちょっとこれを次のように整形してみよう。

　 $A_{000}$ , $A_{001}$ , $A_{002}$ ,
　 $A_{011}$ , $A_{012}$ ,
　 $A_{022}$ ,
　 $A_{111}$ , $A_{112}$ ,
　 $A_{122}$ ,
　 $A_{222}$

ここで、 $i = 0$ の、 $j$ , $k$ 部分に注目すると、添字が２個の場合の上三角行列の添字が並んでいる事に気づくだろう。同様に $i = 1$ の場合は、「1 オリジンの場合の」添字が２個の場合の上三角行列の添字が並んでいる。言い換えると、 $j-i$ , $k-i$ を考えれば、それが $N = 1$ , $M = 2$ の場合の上三角行列の添字になっているということである。このことを考えれば、 $i$ を固定した場合のとりうる添字の数は、 $N^\prime = N-i$ , $M^\prime = M-1$ の要素数の総数に等しいことになる。従って、

　 $U(i, i, ..., i: N, M)$
　　 $= \sum_{n=0}^{i-1} [U(N-n, N-n, ..., N-n: N-n, M-1)+1]$
　 $U(i, j, ..., k: N, M)$
　　 $= U(i, i, ..., i: N, M) + U(j-i, ..., k-i: N-i, M-1)$

となる。即ち、添字の個数が $M-1$ の式さえ分かっていれば、添字の個数が $M$ の式が出せるわけである。例えば $M = 2$ の時は

　 $U(i, j: N, 2)$
　　 $= \sum_{n=0}^{i-1} [U(N-n: 1)+1] + U(j-i: N-i, 1)$
　　 $= \sum_{n=0}^{i-1} (N+1-n) + j-i$
　　 $= i(N+1) - \frac{i(i-1)}{2} + j-i$
　　 $= \frac{i(2N+1-i)}{2} + j$

と求まる。 $I = N-i$ , $J = N-j$ と置き換えると、

　 $U(I, J: N, 2)$
　　 $= \frac{N(N+1)}{2} - \frac{I(I+1)}{2} + N - J$
　　 $= \left(\begin{matrix}N+1 \cr 2\end{matrix}\right) + \left(\begin{matrix}N \cr 1\end{matrix}\right) - \left(\begin{matrix}I+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}J \cr 1\end{matrix}\right)$
　　 $= \left(\begin{matrix}N+1 \cr 2\end{matrix}\right) + \left(\begin{matrix}N \cr 1\end{matrix}\right) + \left(\begin{matrix}N \cr 0\end{matrix}\right) - 1 - \left(\begin{matrix}I+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}J \cr 1\end{matrix}\right)$
　　　(∵ $\left(\begin{matrix}N \cr 0\end{matrix}\right) = 1$ )
　　 $= \left(\begin{matrix}N+1 \cr 2\end{matrix}\right) + \left(\begin{matrix}N+1 \cr 1\end{matrix}\right) - 1 - \left(\begin{matrix}I+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}J \cr 1\end{matrix}\right)$
　　　(∵ パスカルの三角形)
　　 $= \left(\begin{matrix}N+2 \cr 2\end{matrix}\right) - 1 - \left(\begin{matrix}I+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}J \cr 1\end{matrix}\right)$

となる（ $\left(\begin{matrix}a \cr b\end{matrix}\right)$ は $a$ 個の中から $b$ 個を選び取る組み合わせの数（コンビネーション））。M = 3 の場合は、具体的な計算式は省略するが、

　 $U(i, j, k: N, 3)$
　　 $= \frac{i[i^2-3(N+1)i+3N(N+2)+2]}{6} + \frac{j(2N+1-j)}{2} + k$

となる。 $I = N-i$ , $J = N-j$ , $K = N-k$ と置き換えると、

　 $U(I, J, K: N, 3)$
　　 $= \frac{N(N+1)(N+2)}{6} - \frac{I(I+1)(I+2)}{6}$
　　　 $+ \frac{N(N+1)}{2} - \frac{J(J+1)}{2} + N - K$
　　 $= \left(\begin{matrix}N+2 \cr 3\end{matrix}\right) + \left(\begin{matrix}N+1 \cr C 2\end{matrix}\right) + \left(\begin{matrix}N \cr 1\end{matrix}\right)$
　　　 $- \left(\begin{matrix}I+2 \cr 3\end{matrix}\right) - \left(\begin{matrix}J+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}K \cr 1\end{matrix}\right)$
　　 $= \left(\begin{matrix}N+3 \cr 3\end{matrix}\right) - 1 - \left(\begin{matrix}I+2 \cr 3\end{matrix}\right) - \left(\begin{matrix}J+1 \cr 2\end{matrix}\right) - \left(\begin{matrix}K \cr 1\end{matrix}\right)$

となる。従って、一般形は

　 $U(I, J, ..., K: N, M)$
　　 $= \left(\begin{matrix}N+M \cr M\end{matrix}\right) - 1$
　　　 $- \left(\begin{matrix}I+M-1 \cr M\end{matrix}\right) - \left(\begin{matrix}J+M-2 \cr M-1\end{matrix}\right) - ... - \left(\begin{matrix}K \cr 1\end{matrix}\right)$

即ち、

　 $U(i, j, ..., k: N, M)$
　　 $= \left(\begin{matrix}N+M \cr M\end{matrix}\right) - 1$
　　　 $- \left(\begin{matrix}N+M-1-i \cr M\end{matrix}\right) - \left(\begin{matrix}N+M-2-j \cr M-1\end{matrix}\right) - ... - \left(\begin{matrix}N-k \cr 1\end{matrix}\right)$

となる。このことは、

　 $U(i: N, 1)$
　　 $= \left(\begin{matrix}N+1 \cr 1\end{matrix}\right) - 1 - \left(\begin{matrix}N-i \cr 1\end{matrix}\right)$
　　 $= i$
　 $U(i, i, ..., i: N, M+1)$
　　 $= \sum_{n=0}^{i-1} [U(N-n, N-n, ..., N-n: N-n, M)+1]$
　　 $= \sum_{n=0}^{i-1} \left(\begin{matrix}N+M-n \cr M\end{matrix}\right)$
　　　( $\left(\begin{matrix}X \cr X+1\end{matrix}\right) = 0$ と定義する。これは
　　　　　 $\left(\begin{matrix}X \cr X+1\end{matrix}\right) = \frac{X(X-1)...(X-X)}{(X+1)!} = 0$
　　　と矛盾しない)
　　 $= \left(\begin{matrix}N+M \cr M\end{matrix}\right) + \left(\begin{matrix}N+M-1 \cr M\end{matrix}\right) + ... + \left(\begin{matrix}N+M+1-i \cr M\end{matrix}\right)$
　　 $= \left(\begin{matrix}N+M \cr M\end{matrix}\right) + \left(\begin{matrix}N+M-1 \cr M\end{matrix}\right) + ... + \left(\begin{matrix}N+M+1-i \cr M\end{matrix}\right)$
　　　 $+ \left(\begin{matrix}N+M+1-i \cr M+1\end{matrix}\right) - \left(\begin{matrix}N+M+1-i \cr M+1\end{matrix}\right)$
　　 $= \left(\begin{matrix}N+M+1 \cr M+1\end{matrix}\right) - \left(\begin{matrix}N+M+1-i \cr M+1\end{matrix}\right)$
　 $U(j-i, ..., k-i: N-i, M)$
　　 $= \left(\begin{matrix}N+M-i \cr M\end{matrix}\right) - 1 - \left(\begin{matrix}N+M-1-j \cr M\end{matrix}\right) - ... - \left(\begin{matrix}N-k \cr 1\end{matrix}\right)$
　 $U(i, j, ..., k: N, M+1)$
　　 $= \left(\begin{matrix}N+M+1 \cr M+1\end{matrix}\right) - \left(\begin{matrix}N+M+1-i \cr M+1\end{matrix}\right)$
　　　 $+ \left(\begin{matrix}N+M-i \cr M\end{matrix}\right) - 1 - \left(\begin{matrix}N+M-1-j \cr M\end{matrix}\right) - ... - \left(\begin{matrix}N-k \cr 1\end{matrix}\right)$
　　 $= \left(\begin{matrix}N+M+1 \cr M+1\end{matrix}\right) - 1$
　　　 $- \left(\begin{matrix}N+M-i \cr M+1\end{matrix}\right) - \left(\begin{matrix}N+M-1-j \cr M\end{matrix}\right) - ... - \left(\begin{matrix}N-k \cr 1\end{matrix}\right)$
　　　(∵ $\left(\begin{matrix}N+M-i \cr M+1\end{matrix}\right) + \left(\begin{matrix}N+M-i \cr M\end{matrix}\right) = \left(\begin{matrix}N+M+1-i \cr M+1\end{matrix}\right)$ )

のように、数学的帰納法で証明できる。
　超下三角行列の場合も状況は似ている。N = 2, M = 3 の場合、添字は

　 $A_{000}$ , $A_{100}$ , $A_{110}$ , $A_{111}$ , $A_{200}$ , $A_{210}$ , $A_{211}$ , $A_{220}$ , $A_{221}$ , $A_{222}$

と並ぶが、これを

　 $A_{000}$ ,
　 $A_{100}$ ,
　 $A_{110}$ , $A_{111}$ ,
　 $A_{200}$ ,
　 $A_{210}$ , $A_{211}$ ,
　 $A_{220}$ , $A_{221}$ , $A_{222}$

と整形すれば、 $i$ を固定すれば、 $j$ , $k$ は $N^\prime = i$ , $M^\prime = M-1$ の場合の下三角行列の添字になっていることが分かる。従って、

　 $D(i, 0, ..., 0: M)$
　　 $= \sum_{n=0}^{i-1} [D(n, n, ..., n: M-1)+1]$
　　 $= \sum_{n=0}^{i-1} D(n+1, 0, ..., 0: M-1)$
　 $D(i, j, ..., k: M)$
　　 $= D(i, 0, ..., 0: M) + D(j, ..., k: M-1)$

となる。この級数和は単純な形で求められ、

　 $D(i, 0, ..., 0: M)$
　　 $= \sum_{n=0}^{i-1} D(n+1, 0, ..., 0: M-1)$
　　 $= \left(\begin{matrix}M-1+i \cr M\end{matrix}\right)$

となる。これも、

　 $D(i: 1)$
　　 $= i = \left(\begin{matrix}i \cr 1\end{matrix}\right)$
　 $D(i, 0, ..., 0: M+1)$
　　 $= \sum_{n=0}^{i-1} D(n+1, 0, ..., 0: M)$
　　 $= \sum_{n=0}^{i-1} \left(\begin{matrix}M+n \cr M\end{matrix}\right)$
　　 $= \left(\begin{matrix}M \cr M\end{matrix}\right) + \left(\begin{matrix}M+1 \cr M\end{matrix}\right) + ... + \left(\begin{matrix}M+i-1 \cr M\end{matrix}\right)$
　　 $= \left(\begin{matrix}M+1 \cr M+1\end{matrix}\right) + \left(\begin{matrix}M+1 \cr M\end{matrix}\right) + ... + \left(\begin{matrix}M+i-1 \cr M\end{matrix}\right)$
　　　　(∵ $\left(\begin{matrix}M \cr M\end{matrix}\right) = 1 = \left(\begin{matrix}M+1 \cr M+1\end{matrix}\right)$ )
　　 $= \left(\begin{matrix}M+i \cr M+1\end{matrix}\right)$

のように、数学的帰納法で簡単に証明できる。従って、

　 $D(i, j, ..., k: M)$
　　 $= \left(\begin{matrix}M-1+i \cr M\end{matrix}\right) + \left(\begin{matrix}M-2+j \cr M-1\end{matrix}\right) + ... + \left(\begin{matrix}k \cr 1\end{matrix}\right)$

のようになる。
　 $U$ と $D$ の間には、

　 $U(i, j, ..., k: N, M)$
　　 $= \left(\begin{matrix}N+M \cr M\end{matrix}\right) - 1 - D(N-i, N-j, ..., N-k: M)$

という関係がなりたつ。この式に現れる $\left(\begin{matrix}N+M \cr M\end{matrix}\right)$ はまさしく全要素数であり、 $\left(\begin{matrix}N+M \cr M\end{matrix}\right) - 1$ は最後のインデックスである。また、この式の示す通り、基本的に $D$ より $U$ の方が式は面倒になる（例外は $M = 1$ の時で、この場合 $U(i:N, 1) = D(i:1) = i$ ）。

　最後に、競馬などでよく見られるように、同じ値が二度出てこない場合の三角行列を考える。これは非常に単純で、添字の原点をずらせばいいだけの話である。例えば $N$ 頭立ての馬連やワイドの場合に、オッズ表に並んでる通り

01-02 02-03 03-04 ... (N-1)-N
01-03 02-04 03-05 ...
01-04 02-05 03-06
  :
01-N

と並べて、下向き優先のインデクシングをすることを考える。ここで、１頭目のインデックス $x$ から 1 を引き、２頭目のインデックス $y$ から 2 を引くと、

00-00 01-01 02-02 ... (N-2)-(N-2)
00-01 01-02 02-03 ...
00-02 01-03 02-04
  :
00-(N-2)

のようになる。これを下向き優先でインデクシングするということは、上三角行列のインデクシングをすればいいという事に他ならない。従って、インデックスは $U(x-1, y-2: N-2, 2)$ で求める事ができる。また、ちょっと整形して

01-02   ↓    ↓
01-03 02-03   ↓　　一番下を揃えるように下げてやる
01-04 02-04 03-04
  :
01-N  02-N  03-N　...　(N-1)-N

のようにし、右向き優先でインデクシングすることを考えてみよう。この場合も $x$ から 1 を、 $y$ から 2 を引いて

00-00
00-01 01-01
00-02 01-02 02-02
  :
00-(N-2) 01-(N-2) 02-(N-2) ... (N-2)-(N-2)

としてみれば、下三角行列のインデクシングを利用して、 $D(y-2, x-1: 2)$ で求まる事が容易に分かるだろう。三連複の場合も同様に考える事ができ、 $U(x-1, y-2, z-3: N-3, 3)$ や $D(z-3, y-2, x-1: 3)$ でインデクシングすることができる。