一様にランダムな直交行列の生成

問題設定

直交行列 (回転行列) を一様にランダムに生成したい.

直交行列は画像や姿勢データの距離を保存する変換に利用され, ランダムに生成できるとデータ増強等に役に立つことがある.
しかしクオータニオンは次元が増えると使いにくいし, オイラー角は対称性を持たず一様性を得るのに補正が必要なので使いたくない.
任意の次元で統一された手続きでランダムな直交行列が生成できる方法が欲しい.

結論

(1) 正規分布に従う $d \times d$ 個の乱数を生成して正方行列を得る,
(2) (1) で得た正方行列を QR 分解する,
(3) (2) で得た直交行列を必要に応じて三角行列の対角成分によって補正する,

で標準的な分布について一様にランダムな直交行列が得られる.

QR 分解にグラム・シュミットの直交化を用いれば直交行列の計算だけで済んで三角行列の保持は不要になり, (3) の手続きも省略できる.
ただし直交行列の直交性が不安定になる問題がある.

反転を伴わない回転行列を生成する場合は特殊直交行列と呼ばれる行列に変換する必要があり, 上の手続きで生成した行列の行列式を使ってもう一手間加える.

前提知識

• 問題設定と結論の理解のため
  • 線型代数
  • 確率論
• 以降の解説の理解のため
  • 位相空間論 (ハウスドルフ空間, 局所コンパクトぐらいまで)
  • 測度論 (ラドン測度ぐらいまで)
  • 表現論 (位相群, ハール測度の周辺の話題)

位相空間に不慣れなら出て来る空間は全てユークリッド空間だと思って, 測度論が分からなければ「測度 = 確率分布」だと思いながら流し読みしてもいい.
ハール測度は「そういう都合のいい測度 (確率分布) があるんだな」程度に受け取っておこう.

解説

$d$ 次直交行列の集合を直交群と呼び, $O(d)$ と書く.
大筋は次の論文の要約になるが少し流れが違う.

cf. Francesco Mezzadri, "How to generate random matrices from the classical compact groups", 27 Feb 2007, arXiv:math-ph/0609050v2.

一般 (化し過ぎた) 論

標準的な一様分布は $d$ 次元ユークリッド空間の有界な部分集合 $[0, 1]^d$ 上の確率測度であり, ボレル集合族 $\mathcal{B}\left([0, 1]^d\right)$ に制限したルベーグ測度 $\lambda$ を使って,

$$ P(x \in X) := \lambda(X), \quad X \in \mathcal{B}\left([0, 1]^d\right), $$

と書いて定義される.
これを確率空間 $(\Omega, \mathcal{F}, \mu)$ 上に一般化すると, $\Omega$ 上の $\mu$ に関する一様分布を,

$$ P(x \in X) := \mu(X), \quad X \in \mathcal{F}, $$

と定義することができる.
つまり適当に乱数生成した直交行列は一様にランダムな直交行列と見なせる.

で終わってもつまらないので直交群上の代表的な測度であるハール測度を取り上げて話を進める.

ハール測度

以下, 群の作用は左作用であるとする.

局所コンパクト空間 $X$ への局所コンパクト群 $G$ の連続な作用 $G \curvearrowright X$ を与えた時, $B \in \mathcal{B}(X), g \in G$ に対して $gB \in \mathcal{B}(X)$ となる.
この時, $X$ のラドン測度 $\mu$ で任意の $B \in \mathcal{B}(X), g \in G$ に対して,

$$ \mu(gB) = \mu(B), $$

となるものを (左) $G$ 不変であるという.
局所コンパクト群 $G$ には自身への自然な作用 $G \curvearrowright G$ について不変なラドン測度が一意に存在することが知られており, これを $G$ 上の (左) ハール測度と呼ぶ.

$G$ がコンパクト群の場合, ハール測度 $\mu$ は有限測度になることが知られており,

$$ \mu'(B) := \frac{\mu(B)}{\mu(G)}, $$

によって正規化することができてこれは $G$ 上の確率測度となる.

直交群 $O(d)$ はコンパクト群なので有限なハール測度が存在し, これを正規化して一様分布だと見なすことができる.
ハール測度がどういう意味で一様なのかについては議論の余地があり, 論文では球面上一様分布に分解できることを一様性の解釈としている¹.

以下, 有限ハール測度は正規化されたものを考えて, $O(d)$ 上のハール測度に従うランダムな直交行列を生成する方法について述べていこう.

同変写像による乱数の変換

一般にコンパクト群 $G$ 上の確率測度に従う乱数を直接的に生成することは難しいため, 乱数生成可能な空間 $X$ で得た乱数を変換することを考える.

写像による乱数と測度の変換

まず, $X$ 上の測度 $\mu$ から $G$ 上の測度を構成する方法を紹介する.
可測空間 $(X, \mathcal{F}_X), (Y, \mathcal{F}_Y)$ と $X$ 上の測度 $\mu$, 可測写像 $\varphi: X \rightarrow Y$ に対して,

$$ \mu_*\varphi: \mathcal{F}_Y \ni F \mapsto \mu(\varphi^{-1}(F)) \in \mathbb{R}, $$

は $(Y, \mathcal{F}_Y)$ 上の測度になり, $\varphi$ による $\mu$ の像測度と呼ばれる.
定義から, $\mu$ が $X$ 上の確率測度ならば像測度 $\mu _*\varphi$ は $Y$ 上の確率測度となることに注意.

像測度について, $Y$ の可積分関数 $f: Y \rightarrow \mathbb{R}$ に対して,

$$ \int_Xf\circ\varphi d\mu = \int_Yfd\mu_*\varphi, $$

となることが知られている.
これは, $x_0, x_1, \dots, x_{n-1} \in X$ が確率測度 $\mu$ に従う点列ならば, $\varphi(x_0), \varphi(x_1), \dots, \varphi(x_{n-1}) \in Y$ が $\mu_*\varphi$ に従う点列であることを表している.
つまり, $X$ の可積分関数 $f'$ の $\mu$ に関する積分 (期待値) を,

$$ \frac{1}{n}\sum_{i=0}^{n-1}f'(x_i) \sim \int_Xf'd\mu, $$

と近似できるならば, $Y$ の可積分関数 $f$ の $\mu_*\varphi$ に関する積分について,

$$ \frac{1}{n}\sum_{i=0}^{n-1}f(\varphi(x_i)) \sim \int_Xf\circ\varphi d\mu = \int_Yfd\mu_*\varphi, $$

という近似が成り立つことが示される².

不変なラドン測度の構成

続いて $X$ から $G$ への可測写像 $\varphi$ に対して $\mu_*\varphi$ が $G$ 上のハール測度となる条件を確認しよう.
$G$ 上のハール測度とは $G$ 上の $G$ 不変なラドン測度のことであった.

下のサイト等により, 第二可算な局所コンパクトハウスドルフ空間上の (ボレル集合族で定義される) 有限測度はラドン測度となることが分かっている.
よって, $Y$ が第二可算な局所コンパクトハウスドルフ空間の時, 可測写像 $\varphi: X \rightarrow Y$ による $X$ 上の確率測度の像測度は $Y$ 上の確率測度で有限であり, ラドン測度となることが分かる.
cf. 測度と積分7：局所コンパクトHausdorff空間上のRadon測度 - Mathpedia

$G$ 不変な像測度を得るためには, 次で定義される同変写像を利用するといい.
局所コンパクト群 $G$ が作用する 2 つの空間 $X, Y$ 間の連続写像 $\varphi: X \rightarrow Y$ で,

$$ \varphi(gx) = g\varphi(x), \quad x \in X, \ g \in G, $$

となるものを (左) $G$ 同変であるという.
この時, $S \subset X, g \in G$ に対して $\varphi^{-1}(gS) = g\varphi^{-1}(S)$ に注意して, $\mu$ を $X$ 上の $G$ 不変なラドン測度とすると, $B \in \mathcal{B}(Y)$ に対して,

$$ \mu_* \varphi(gB) = \mu(g\varphi^{-1}(B)) = \mu(\varphi^{-1}(B)) = \mu_*\varphi(B), $$

となるため $\mu_*\varphi$ は $Y$ 上の $G$ 不変なラドン測度となる.

以上の議論で $Y = G$ と置けば, $X, G$ が共に第二可算な局所コンパクトハウスドルフ空間の時, $X$ の $G$ 不変な確率測度 $\mu$ の $G$ 同変写像 $\varphi$ による像測度は $G$ 上のハール測度となることが導かれる.

以上のことから, 次のように直交群 $O(d)$ 上のハール測度に従う乱数を生成できることが分かった.
つまり, $O(d)$ は $n$ 次元ユークリッド空間 $\mathbb{R}^n$ に行列のベクトルへの積によって作用し, これらの集合は第二可算な局所コンパクトハウスドルフ空間なので, $\mathbb{R}^n$ 上の $O(d)$ 不変な有限ラドン測度に従う乱数を $O(d)$ 同変な写像で変換すると $O(d)$ 上のハール測度に従う乱数が得られる.
$O(d)$ と $\mathbb{R}^n$ で次元の記号を分けているのは $n = d \geq 2$ の時 $\mathbb{R}^d$ から $O(d)$ への $O(d)$ 同変な写像が存在しないため (演習).

一般線型群から直交群への直交群の作用について同変な写像

$d$ 次直交群 $O(d)$ は, 正方行列のうち各列が $\mathbb{R}^d$ の正規直交基底を成すもの全体として与えられるので, 正規直交基底の生成方法から $\mathbb{R}^{d\times d} \rightarrow O(d)$ の変換方法を考えてみる.

有名な方法として, 正則行列から次の手続きによって正規直行基底を得るグラム・シュミット正規直交化 $\Gamma: GL(d) \rightarrow O(d)$ がある.
正則行列の集合である $GL(d)$ (一般線型群と呼ばれる) は行列式が $0$ の行列を含まないので $\mathbb{R}^{d\times d}$ とは異なる空間じゃないの, という点については後述.

$A \in GL(d)$ の各列が $a_0, a_1, \dots, a_{d-1} \in \mathbb{R}^d$ からなるとして ($A = [a_0, a_1, \dots, a_{d-1}]$ と書く), 以下の操作を行う.
(0) $i := 0$,
(1) $\sigma_i(A) := a_i-\sum _{j=0}^{i-1}\langle a_i \vert \gamma_j(A) \rangle \gamma_j(A)$ ($i = 0$ なら右辺第二項は $0$),
(2) $\gamma_i(A) := \sigma_i(A)/\|\sigma_i(A)\|$,
(3) $i=d-1$ なら終了, そうでないなら $i \leftarrow i+1$ して (1) へ.
こうして $\Gamma(A) := [\gamma_0(A), \gamma_1(A), \dots, \gamma _{d-1}(A)]$ を得る操作をグラム・シュミット直交化と呼び, $\Gamma(A) \in O(d)$ となる.

$\Gamma$ は連続写像の合成で書けるので連続で, 直交群の左作用について同変になる³.
証明は省くが, $P \in O(d)$ を取って $A \leftarrow PA$ を代入して, 直交行列が内積とノルムを保存することを利用して式変形していくと示せる.
よって, $GL(d)$ 上 $O(d)$ 不変な確率測度 $\mu$ とそれに従う乱数が得られれば, それを $\Gamma$ で写したものは $O(d)$ 上のハール測度に従う乱数となる.

ただし, グラム・シュミット直交化は数値計算の安定性に問題があり, 論文ではハウスホルダー変換による直交化と三角行列による補正によってランダムな直交行列を生成する方法を採用している.
この記事では理論的なところを重視しているので最後に軽く触れる程度に留める.

$GL(d)$ 上 $O(d)$ 不変な有限ラドン測度

$GL(d)$ 上 $O(d)$ 不変な有限ラドン測度の構成には $O(d) \curvearrowright \mathbb{R}^d$ の等長性を利用すると良い.
$GL(d)$ 上の確率測度で確率密度関数が $p(a_0, a_1, \dots, a_{d-1}) = \hat{p}(\|a_0\|, \|a_1\|, \dots, \|a_{d-1}\|)$ と書けるものを考えると, $P \in O(d)$ に対して,

$$ \begin{align}
p(P[a_0, a_1, \dots, a_{d-1}]) & = p(Pa_0, Pa_1, \dots, Pa_{d-1}) \\
& = \hat{p}(\|Pa_0\|, \|Pa_1\|, \dots, \|Pa_{d-1}\|) \\
& = \hat{p}(\|a_0\|, \|a_1\|, \dots, \|a_{d-1}\|) \\
& = p(a_0, a_1, \dots, a_{d-1}),
\end{align} $$

より $p$ に対応する確率測度は $O(d)$ 不変になる.
($\mathbb{R}^d$ 上の) 確率密度関数と確率測度の $O(d)$ 不変性が同値になるのはルベーグ測度の変数変換公式等による.

この条件を満たす確率密度関数を持つ $GL(d)$ 上の確率測度としては, $d^2$ 次元標準正規分布が挙げられる.
正規分布はそれに従う乱数の生成方法がいくつか知られていて多くの乱数ライブラリにも実装されているので, 変換元となる確率測度として都合が良い.
一般的な方法だと正方行列空間 $M(d) = \mathbb{R}^{d \times d}$ 上で乱数を生成してしまう点が気になるところだが, $GL(d)$ は $M(d)$ で稠密で $M(d) \setminus GL(d)$ の測度は $0$ になるので無視することができ, 実際, 以下のコードは実質的に無限ループになる.
例外が心配ならコード上で $\det{A} = 0$ の時に再生成するコードを差し込んでおけば良い.

import numpy as np
rng = np.random.default_rng()

d = 2
while True:
    A = rng.normal(size=[d, d])
    if np.linalg.det(A)==0:
        break

$O(d)$ 上ハール測度に関する一様分布に従う乱数の生成

以上の議論から,
(1) $d^2$ 次元正規分布に従う乱数を生成,
(2) 生成した乱数をグラム・シュミット正規直交化,
の手続きで直交群上の標準的な測度に関する一様にランダムな直交行列を生成できることが分かった.

Python で NumPy を使う場合, QR 分解を実行する numpy.linalg.qr では論文で扱われているハウスホルダー変換による直交化が行われてしまうため, グラム・シュミット正規直交化を自前で実装しよう.
numpy.linalg.qr では 2 次元配列しか渡せないため自前実装では複数のランダムな直交行列を生成するのが捗るという利点もある.
生成された乱数が直交群上のハール測度に従っていることの確認は難しいので略.

import numpy as np
rng = np.random.default_rng()

# 正規分布に従う GL(d) の乱数を n 個生成
def normalGL(n, d):
    # det(A)==0 の行列が生成されなくなるまで繰り返し
    while True:
        A = rng.normal(size=[n, d, d])
        crit = np.linalg.det(A)!=0
        if np.all(crit): break
    return A

# グラム・シュミット正規直交化
def GramSchmidt(A):
    k = np.shape(A)[-1]
    Gmm = np.array(A)
    for i in range(k):
        ip = np.sum(Gmm[..., i:i+1]*Gmm[..., :i], axis=-2, keepdims=True)
        sgm = Gmm[..., i]-np.sum(ip*Gmm[..., :i], axis=-1)
        Gmm[..., i] = sgm/np.linalg.norm(sgm, axis=-1, keepdims=True)
    return Gmm

# 一様にランダムな直交行列の生成
def uniformO(n, d):
    return GramSchmidt(normalGL(n, d))

# 10 次直交行列を 10**4 個生成
n = 10**4
d = 10
P = uniformO(n, d)

付録 1. 回転行列を生成する場合

画像の回転や姿勢制御に用いられる回転行列を生成する場合, 直交行列 $P \in O(d)$ のうち行列式が $\det{P} = 1$ となるものを取る必要がある.
これを特殊直交行列と呼び, その集合を $SO(d)$ と書いて特殊直交群と呼ぶ.

$SO(d)$ もコンパクト群なので有限ハール測度を持ち, $O(d)$ 上の $SO(d)$ 不変測度と $O(d)$ から $SO(d)$ への $SO(d)$ 同変写像から $SO(d)$ 上のハール測度に従う乱数を生成できる.

不変測度については, $SO(d) \subset O(d)$ から $O(d)$ 不変測度が $SO(d)$ 不変測度にもなるため, $O(d)$ 上のハール測度を選ぶことができる.
同変写像 $\varphi: O(d) \rightarrow SO(d)$ は以下で構成される.

$$ \begin{align}
\Delta(P) & := \mathrm{diag}(\det{P}, 1, \dots, 1) \\
& = \begin{bmatrix}
\det{P} & 0 & \dots & 0 \\
0 & 1 & \ddots & \vdots \\
\vdots & \ddots & \ddots & \vdots \\
0 & \dots & 0 & 1 \\
\end{bmatrix}, \\
\varphi(P) & := P\Delta(P). \\
\end{align} $$

$\mathrm{diag}(\lambda_0, \lambda_1, \dots \lambda_{d-1})$ は対角成分に $\lambda_0, \lambda_1, \dots \lambda_{d-1}$ を持つ対角行列である.

これにより, 上の手続きで生成したランダムな直交行列 $P$ に $\Delta(P)$ を右から掛ければランダムな回転行列となることが分かる⁴.

# 一様にランダムな特殊直交行列 (回転行列) の生成
# 直交行列 P を生成してから Delta(P) を掛ける
def uniformSO(n, d):
    P = uniformO(n, d)
    D = np.ones_like(P)*np.eye(d)
    D[..., 0, 0] = np.linalg.det(P)
    return P@D

# 10 次特殊直交行列を 10**4 個生成
n = 10**4
d = 10
P = uniformSO(n, d)

付録 2. ハウスホルダー変換その他の直交化法を使った場合

詳細は論文参照.
$d^2$ 次元正規分布に従って生成した $A \in GL(d)$ を $A = QR$ と QR 分解して三角行列成分 $R$ の対角成分 $r_0, r_1, \dots, r_{d-1}$ を取り出して,

$$ D := \mathrm{diag}\left(\frac{r_0}{|r_0|}, \frac{r_1}{|r_1|}, \dots, \frac{r_{d-1}}{|r_{d-1}|}\right), $$

と置く.
すると $P = QD$ と補正すれば $P$ はグラム・シュミット正規直交化 $\Gamma(A)$ に一致し, 従って $O(d)$ 上のハール測度に従う乱数になる.
$A \in GL(d)$ を仮定しているので $R$ の対角成分は非ゼロとなっている.

これによりハウスホルダー変換による方法を初めとした, グラム・シュミット正規直交化より安定的な手法を使ってランダムな直交行列が得られる.

1. 多様体 (計量) に疎いので確認できていないが直交群のハール測度はリーマン測度の定数倍に一致する気がする. $O(d)$ の $O(d)$ への作用は $O(d \times d)$ で書けることから $O(d)$ 上のリーマン測度の $O(d)$ 不変性が示せないだろうか. ↩

2. 話が飛ぶようだが, この近似は点列が確率測度に従うことの条件として扱われる. 大数の法則において可測関数 (確率変数) の列は確率測度に従う点列で置き換えられるという解釈の元でそのように特徴付けられる. ↩

3. 右作用については同変にならない. これは QR 分解 $A = \Gamma(A)R$ を考えた時, $A$ への右作用が直接的には三角行列成分 $R$ に作用することに対応すると考えられる. ↩

4. 実は左から掛けてもいい. $\Delta(P)$ を左から掛ける写像は左同変でなく右同変となるが, ハール測度は右不変でもあるので結局はハール測度に従うランダムな行列が得られる. ↩