行列の成分座標と配列のインデックス番号の関係

行列やグラフを扱うプログラムを書いていると,M x N 行列の成分座標を表す整数の組 (i, j) と,要素の通し番号 t の変換をする場合がしばしばある.たとえば,C プログラム中で行列を1次元配列として扱う場合,(i, j) 成分にアクセスするためにはこれをインデックス番号に変換する必要があるし,その逆もときどき必要になったりする.

M x N 行列のすべての成分を1次元配列に収める場合は簡単である.成分座標 (i, j) からインデックス番号 t,あるいはその逆変換は,次のようにして行える.ここでは行列の配置は行ベースで,成分座標もインデックス番号も 0-origin とする.

t <- i + N * (j - 1)        # (i, j) -> t

i <- t %% N                 # t -> (i, j)
j <- floor(t / N)

ところで,N x N 対称行列の場合は,行列の右上側だけを保持することで,記憶容量を節約したり,プログラムミスによる対称条件の破壊を防ぐことができる((i, j) 成分を更新したのに (j, i) 成分を更新するのを忘れる,など).このために,i <= j を満たすような (i, j) 成分の要素を辞書順で配列に格納する,というコードをよく書く.たとえば次のようなコード. [code language="r"] vec <- numeric(length=N * (N + 1) / 2) t <- 0 for (i in 1 : (N - 1)) for (j in i : N) { vec[t] <- calculate_value(i, j) t <- t + 1 } [/code] この場合,要素のアクセスは for や while などを使って逐次的に行うことが多いが,解析的に一発で座標を求める方法はないものか,というのが今回の記事のテーマ. (i, j) から t への変換は,対角成分が三角数になることを利用すればすぐに思いつく.下三角の場合はもっと複雑. [code language="r"] t <- i + j * (j - 1) / 2 # (*) [/code] 問題は t から (i, j) に直すとき.1変数から2変数への変換なので,単純に式変形,というわけにはいかない.結論から言うと,次のようになる. [code language="r"] i <- t - j * (j - 1) / 2 j <- floor((1 + sqrt(8 * t - 7)) / 2) [/code] j が求まれば (*) 式から i を求めるのは簡単なので,問題は j をどうやって求めるかになる.(*) の式を j について解くと,これは2次方程式なので平方根を含んだ解が出てくる.2つの解のうちどちらが真の解かは,j >= 1 の条件からすぐにわかる.このままでは i を含んでいるので,i を含まない形に変形したいのだが,問題の条件を考えれば i を動かしても j の整数部分は変わらないはずなので,i = 1 で固定することで式中から i を消す.固定する i はうまく選ばないと平方根の中が負になったり,j が N より大きくなったりするので注意.また,行ベースや 0-origin,対角成分を含まない場合などを考えるときは,式や i の固定位置が少しずつ変わるのでこれも注意.

解に平方根が含まれるのはすぐ予想できるけど,最後の部分の厳密な証明はぱっと思いつかなかった.ほかにこんなこと考えてる人いないかしら?

2013-07-12 修正: 行ベース・0 – origin の場合を書いていたけど,R のコードを書いてるなら列ベース・ 1-origin のほうがいいんじゃないか?と思ったので修正.

コメントを残す

メールアドレスが公開されることはありません。