問題設定

$d$ 次元空間上でベクトルの長さを変えないよう変換する行列を直交行列と呼ぶ.
$3$ 次元空間で回転行列と呼ばれる「向きを保存する」直交行列を扱う場合, クオータニオン (四元数) を用いて記述する方法があり, 直交行列を回転軸と回転平面に分解する操作に密接に関わっている.

それで, その直交行列の回転軸と回転平面への分解を一般次元で論じられないかという問題.

「回転行列」という用語について

一般の $d$ 次元空間において「回転行列」と言った場合,
(a) 行列式が $1$ である $d$ 次直交行列 (特殊直交行列),
(b) $2$ 次元平面上のベクトルを回転させ, その平面に直交する $d-2$ 次元成分を保存する直交行列の一種,
の二つの意味があり, $d \leq 3$ では両者は一致する.

(a) は (b) を含むためこの記事では (a) の意味で回転行列の一般次元化 (さらに行列式が $-1$ の場合への拡張) とその分解を考えるが, (b) の意味では次の記事やそこで紹介されている論文で一般次元化と分解について述べられている.

結論

(1) $P \in O(d)$ が実固有値を持つ場合

$d$ 次直交行列 $P \in O(d)$ ($O(d)$ は直交行列の集合) が実固有値 $\sigma$ を持つ時, ある $R \in O(d), P' \in O(d-1)$ によって,

$$ P = R \begin{bmatrix} \sigma & 0 \\ 0 & P' \end{bmatrix} R^t, $$

と分解できる.

直交行列の性質から $\sigma = \pm1$ であり, $R$ の第 $1$ 列を $v$ と置くとこれは $\sigma$ に対応する固有ベクトルになる.
この行列の分解は $P$ の作用が $v$ を回転軸として直交する超平面を $P'$ で回転させる操作として表されることを表している.

(2) $P \in O(d)$ が実でない固有値を持つ場合

$d \geq 2$ の時, $d$ 次直交行列 $P \in O(d)$ が (実でない) 複素固有値 $\lambda$ を持つ時, $P$ は,

$$ P = R\begin{bmatrix} \rho(\lambda) & 0 \\ 0 & P' \end{bmatrix}R^t, $$

という形に分解できる¹.
ただし, $R \in O(d), P' \in O(d-2)$ で, $\lambda = \alpha+i\beta$ ($\alpha, \beta \in \mathbb{R}$) と書けるとすると,

$$ \rho(\lambda) = \begin{bmatrix} \alpha & -\beta \\ \beta & \alpha \end{bmatrix}, $$

である.
直交行列の性質から $|\lambda| = 1$ であることに注意すると, $\rho(\lambda)$ は, $2$ 次元ベクトルを $\lambda$ の偏角だけ回転させる行列と言い換えることができる.

この分解は, $P$ による作用が,
(a) $R$ の左 $2$ 列のベクトルが張る平面 (回転平面) について $\lambda$ の偏角だけ回転させる,
(b) $R$ の右 $d-2$ 列のベクトルが張る超平面に $P'$ を作用させる,
という操作に分解されることを表している.
$R$ の左 $2$ 列のベクトルは $\lambda$ に対応する固有ベクトルの虚部と実部からなる.

(3) 一般の直交行列の分解

任意の $d$ 次直交行列 $P \in O(d)$ は (1), (2) を繰り返し適用することで,

$$ P = R \begin{bmatrix}
\sigma_1 & & & & & 0 \\
& \ddots & & & & \\
& & \sigma_{d'} & & & \\
& & & \rho(\lambda_1) & & \\
& & & & \ddots & \\
0 & & & & & \rho(\lambda_{d''})
\end{bmatrix} R^t, $$

と分解できる.
右辺中央の行列の非対角 (非ブロック対角) 成分は全て $0$ で, $\sigma_i$ は $P$ の実固有値, $\lambda_i$ は実でない固有値である.
$R$ の各列 (つまり固有ベクトル) は $P$ による作用の回転軸または回転平面の基底を与える.

前提知識

線型代数, 特に直交行列と固有値周りの話題.
ブロック行列が多用されるので記法に慣れているといい.
クオータニオンを使った回転行列の記述が分かっているとイメージが捗るかもしれない.

解説

証明は端折り気味なので後で補完記事を書くかもしれない.

直交行列

$d$ 次正方行列 $P \in M(d, \mathbb{R})$ で, ($\mathbb{R}^d$ の標準内積に関する) 転置 $P^t$ について $P^{-1} = P^t$ が成り立つものを ($d$ 次の) 直交行列と呼び, $d$ 次直交行列全体の集合を $O(d)$ と書く.

直交行列の性質は次のサイト等を参照.
この記事で重要なものとしては直交行列の各列または各行が正規直交基底をなすことと, 全ての固有値の絶対値が $1$ になることの二つ.

クオータニオンと回転行列の分解

直交行列 $P \in O(d)$ で $\det{P} = 1$ となるものを回転行列 (または特殊直交行列) と呼ぶ.
$d = 3$ の場合, 回転行列はクオータニオン²と呼ばれる長さ $1$ の $4$ 次元ベクトルと相互に変換できることが知られている.
詳細は以下の記事等から.

クオータニオン $q \in \mathbb{R}^4$ は単位ベクトル $v = [v_0, v_1, v_2]^t \in \mathbb{R}^3$ ($\|v\| = 1$) と $\theta \in \mathbb{R}$ によって,

$$ q = \begin{bmatrix} v_0\sin{\frac{\theta}{2}} \\ v_1\sin{\frac{\theta}{2}} \\ v_2\sin{\frac{\theta}{2}} \\ \cos{\frac{\theta}{2}}, \end{bmatrix} $$

と書けて, これは $v$ を回転軸として $v$ に直交する平面上で $\theta$ だけ回転させる直交行列に対応する.
この時, $P$ は第 $1$ 列が $v$ である回転行列 $R$ を用いて,

$$ P = R \begin{bmatrix}
1 & 0 & 0 \\
0 & \cos{\theta} & -\sin{\theta} \\
0 & \sin{\theta} & \cos{\theta} \\
\end{bmatrix} R^t, $$

と分解される.
$R$ の取り方は第 $1$ 列が $v$ であれば任意である.

この分解は,
(1) $R$ の左 $1$ 列に回転軸が表れる分解,
(2) $R$ の右 $2$ 列に回転平面の基底が表れる分解,
という二つの側面があり, これが一般次元の直交行列にどう拡張されるかを見ていく.

一般次元の直交行列の分解

(1) $P \in O(d)$ が実固有値を持つ場合

$\sigma \in \mathbb{R}$ を $P \in O(d)$ の実固有値として, $v \in \mathbb{R}^d$ を対応する単位固有ベクトル ($\|v\| = 1$) とする.
この時, 第 $1$ 列が $v$ である直交行列 $R \in O(d)$ が取れて, $R = \begin{bmatrix} v & Q \end{bmatrix}$ とおくと $R^tR = I_d$ ($I_d$ は $d$ 次単位行列) より $Q^tv = 0$ となる.
また $\sigma$ は $P^t$ の固有値でもあって $v$ が対応する固有ベクトルとなることから,

$$ R^tPR = \begin{bmatrix} v^t \\ Q^t \end{bmatrix} P \begin{bmatrix} v & Q \end{bmatrix} = \begin{bmatrix}
v^tPv & v^tPQ \\
Q^tPv & Q^tPQ \\
\end{bmatrix} = \begin{bmatrix}
\sigma & 0 \\
0 & Q^tPQ \\
\end{bmatrix}, $$

となる.

よって, $R^tPR$ が直交行列で $P' = Q^tPQ$ と置くとこれも直交行列となるので,

$$ P = R \begin{bmatrix}
\sigma & 0 \\
0 & P' \\
\end{bmatrix} R^t, $$

という分解が得られる.

(2) $P \in O(d)$ が実でない固有値を持つ場合

$P \in O(d)$ が実でない固有値を持つ時の実空間での振る舞いについては次の記事を参照.

$\lambda \in \mathbb{C} \setminus \mathbb{R}$ を $P$ の実でない固有値として $v \in \mathbb{C}^d$ を対応する長さ $\sqrt{2}$ の固有ベクトルとする.
この時 $\lambda = \alpha+i\beta$ ($\alpha, \beta \in \mathbb{R}$), $v = t+iu$ ($t, u \in \mathbb{R}^d$) と置くと, 上の記事の結果から,

$$ \begin{cases}
Pt = \alpha t-\beta u, \\
Pu = \beta t+\alpha u, \\
P^tt = \alpha t+\beta u, \\
P^tu = -\beta t+\alpha u, \\
\end{cases} $$

が得られる.

また, $\bar{\lambda} = \alpha-i\beta$ も $P$ の実でない固有値で $\bar{v} = t-iu$ は対応する固有ベクトルとなり, $v$ と $\bar{v}$ は直交する.
これにより $\langle v\vert\bar{v}\rangle = (t+iu)^2 = 0, \|v\|^2 = \|t\|^2+\|u\|^2 = 2$ を解くと $v$ の実部 $t$ と虚部 $u$ はともに長さ $1$ で互いに直交するベクトルとなることが分かる.
よって, 直交行列 $R$ で第 $1$ 列が $u$ で第 $2$ 列が $t$ であるものを取ることができ, $R = \begin{bmatrix} u & t & Q \end{bmatrix}$ と書くと, (1) と同様に,

$$ R^tPR = \begin{bmatrix}
u^tPu & u^tPt & u^tPQ \\
t^tPu & t^tPt & t^tPQ \\
Q^tPu & Q^tPt & Q^tPQ \\
\end{bmatrix}
= \begin{bmatrix}
\alpha & -\beta & 0 \\
\beta & \alpha & 0 \\
0 & 0 & Q^tPQ \\
\end{bmatrix}, $$

が得られて $P' = Q^tPQ$ は直交行列となる.

従って,

$$ P = R \begin{bmatrix}
\rho(\lambda) & 0 \\
0 & P' \\
\end{bmatrix} R^t, $$

と分解できることが分かる.

(3) 一般の直交行列の分解

直交行列 $P \in O(d)$ を (1) によって,

$$ P = R \begin{bmatrix}
\sigma & 0 \\
0 & P' \\
\end{bmatrix} R^t, $$

と分解したとすると, ($P'$ が実固有値を持つと仮定して) $P'$ にさらに (1) を適用すれば,

$$ P = R \begin{bmatrix}
\sigma & 0 \\
0 & R' \begin{bmatrix}
\sigma' & 0 \\
0 & P'' \\
\end{bmatrix} R'^t \\
\end{bmatrix} R^t = R^* \begin{bmatrix}
\sigma & 0 & 0 \\
0 & \sigma' & 0 \\
0 & 0 & P'' \\
\end{bmatrix} R^{*t}, $$

が得られる.
ただし $R^* = R \begin{bmatrix}
1 & 0 \\
0 & R' \\
\end{bmatrix}$ である.
この時 $P'$ の固有値 $\sigma'$ は $P$ の固有値にもなり, これを繰り返すと実固有値を持たない直交行列 $P^*$ によって,

$$ P = R^{**} \begin{bmatrix}
\sigma_1 & & & 0 \\
& \ddots & & \vdots \\
& & \sigma_{d'} & 0 \\
0 & \dots & 0 & P^* \\
\end{bmatrix} R^{**t}, $$

と分解される.
あとは同様に (2) を繰り返し適用すれば,

$$ P = R^{***} \begin{bmatrix}
\sigma_1 & & & & & \\
& \ddots & & & & \\
& & \sigma_{d'} & & & \\
& & & \rho(\lambda_1) & & \\
& & & & \ddots & \\
& & & & & \rho(\lambda_{d''})
\end{bmatrix} R^{***t}, $$

という分解が得られる³.

分解行列 $R \in O(d)$ の構成

与えられた直交する単位ベクトル $v_1, v_2, \dots, v_k$ に対して, $R_{(k)} = \begin{bmatrix} v_1 & v_2 & \dots & v_k & Q_{(k)} \end{bmatrix}$ と書ける直交行列 $R_{(k)} \in O(d)$ を構成する方法を考える ($Q_{(k)}$ は $d \times (d-k)$ の行列).

$R_{(k-1)} = \begin{bmatrix} v_1 & v_2 & \dots & v_{k-1} & Q_{(k-1)} \end{bmatrix}$ と書ける $R_{(k-1)} \in O(d)$ が得られたとすると,

$$ R_{(k-1)}^tR_{(k)} = \begin{bmatrix} v_1^t \\ v_2^t \\ \vdots \\ v_{k-1}^t \\ Q_{(k-1)}^t \end{bmatrix} \begin{bmatrix} v_1 & v_2 & \dots & v_k & Q_{(k)} \end{bmatrix} = \begin{bmatrix}
I_{k-1} & 0 \\
0 & R_{(k)}' \\
\end{bmatrix}, $$

となり (ただし $R_{(k)}' = \begin{bmatrix} Q_{(k-1)}^tv_k & Q_{(k-1)}^tQ_{(k)} \end{bmatrix}$), $R_{(k)}' \in O(d-k+1)$ を構成できれば,

$$ R_{(k)} = R_{(k-1)} \begin{bmatrix}
I_{k-1} & 0 \\
0 & R_{(k)}' \\
\end{bmatrix}, $$

により $R_{(k)}$ が得られることが分かる.
$Q_{(k)}$ には十分な任意性があるため第 $1$ 列が $Q_{(k-1)}^tv_k$ である直交行列の構成法を考えればよく, $k = 1$ の場合に $R_{(1)}$ が得られればいいことが分かる.

$R_{(1)}$ を得る具体的方法としては,
(a) $v_1$ を含む適当な $\mathbb{R}^d$ の基底を取り, $v_1$ を基準に正規直交化する,
(b) $e_1-v_1$ に関するハウスホルダー変換行列を利用する ($e_1$ は第 $1$ 成分のみ $1$ の標準単位ベクトル),
等があり, 初めに挙げた記事で紹介されている $v_1$ の極座標表示から $R_{(1)}$ を生成する方法も利用できる.

実装と動作確認

import numpy as np
np.set_printoptions(precision=4, suppress=True) # 表示桁数を抑える
rng = np.random.default_rng()

直交行列は下の記事に従って生成できる⁴.
奇数次の正方行列は必ず実固有値を持つので確認に便利.

d = 7
Q, R = np.linalg.qr(rng.normal(size=[d, d]))
P = Q@np.diag(np.sign(np.diag(R)))
print(P)
print(np.linalg.eigvals(P)) # 固有値も見ておく

[[ 0.0764 -0.7666  0.1754  0.2687 -0.3725 -0.078   0.3983]
 [ 0.1573 -0.1752 -0.5905 -0.3887  0.3084  0.3093  0.504 ]
 [ 0.0777 -0.1336 -0.7159  0.4976 -0.1575  0.0129 -0.4371]
 [-0.2351  0.0277  0.113  -0.0402 -0.3878  0.8751 -0.1154]
 [ 0.467   0.3444  0.1653  0.6261  0.2291  0.269   0.3452]
 [-0.7849  0.2259 -0.1747  0.3175  0.0083 -0.1202  0.4324]
 [ 0.2719  0.4397 -0.1937 -0.1863 -0.7337 -0.2131  0.2772]]
[ 0.839 +0.5441j  0.839 -0.5441j  0.2673+0.9636j  0.2673-0.9636j
 -0.8407+0.5415j -0.8407-0.5415j -1.    +0.j    ]

初めに (1) の分解の実装と動作確認.
第 $1$ 列が特定の単位ベクトルになる直交行列は np.linalg.qr(v, mode='complete')[0] の一行で生成できるのでそれを使う.

def decomp_1(P):
    l, e = np.linalg.eig(P)
    j = np.where(l.imag==0)[0][0] # 実固有値のインデクスを取得
    v = e[:, j:j+1].real
    R = np.linalg.qr(v, mode='complete')[0]
    D = R.T@P@R
    return R, D

R, D = decomp_1(P)
print(R)
print(D)

[[-0.1378 -0.1928 -0.7629  0.449   0.0624 -0.3952  0.0155]
 [-0.1928  0.9673 -0.1293  0.0761  0.0106 -0.067   0.0026]
 [-0.7629 -0.1293  0.4885  0.301   0.0419 -0.2649  0.0104]
 [ 0.449   0.0761  0.301   0.8228 -0.0246  0.1559 -0.0061]
 [ 0.0624  0.0106  0.0419 -0.0246  0.9966  0.0217 -0.0008]
 [-0.3952 -0.067  -0.2649  0.1559  0.0217  0.8628  0.0054]
 [ 0.0155  0.0026  0.0104 -0.0061 -0.0008  0.0054  0.9998]]
[[-1.     -0.     -0.     -0.      0.      0.      0.    ]
 [ 0.     -0.041  -0.6033 -0.4442  0.3702  0.3313  0.4362]
 [-0.      0.4895 -0.4017  0.0633  0.0569  0.2889 -0.7129]
 [-0.     -0.3069  0.055   0.1407 -0.5244  0.7785  0.0444]
 [ 0.      0.2132 -0.1778  0.8484  0.2376  0.082   0.3742]
 [ 0.      0.6885  0.5413 -0.233   0.0741  0.3094  0.2783]
 [ 0.      0.3808 -0.3834 -0.0696 -0.723  -0.3136  0.2865]]

続いて (2) の分解の実装と動作確認.
R は (1) と同じく np.linalg.qr で生成する.
今度は D の左上 2x2 が $\rho(\lambda)$ の形をしていることを確認しておく.

def decomp_2(P):
    l, e = np.linalg.eig(P)
    j = np.where(l.imag!=0)[0][0] # 非実固有値のインデクスを取得
    v = np.sqrt(2)*e[:, j:j+1]
    t, u = v.real, v.imag
    R = np.linalg.qr(np.hstack([u, t]), mode='complete')[0]
    D = R.T@P@R
    return R, D

R, D = decomp_2(P)
print(R)
print(D)

[[-0.4807 -0.0969  0.0108  0.4433  0.387   0.5945 -0.2444]
 [ 0.3977 -0.4602  0.3213  0.2241 -0.066  -0.1584 -0.6687]
 [-0.0971  0.3179  0.93   -0.0458  0.0171  0.0387  0.1442]
 [ 0.3291  0.3647 -0.0773  0.8222 -0.0948 -0.1324  0.2239]
 [ 0.3743  0.0625 -0.011  -0.1208  0.9012 -0.1503  0.0806]
 [ 0.5941  0.0402 -0.0046 -0.1782 -0.1553  0.7615  0.0987]
 [-0.     -0.7344  0.1602  0.1696  0.0188  0.0006  0.6371]]
[[ 0.839  -0.5441  0.      0.      0.     -0.      0.    ]
 [ 0.5441  0.839  -0.      0.     -0.      0.     -0.    ]
 [ 0.     -0.     -0.923   0.3289 -0.097   0.0801 -0.1552]
 [-0.     -0.     -0.0286 -0.2697 -0.5699  0.7069  0.3193]
 [ 0.     -0.      0.2508  0.8198  0.2142  0.3609  0.2983]
 [-0.      0.     -0.0676  0.1461 -0.4209 -0.5893  0.6706]
 [-0.     -0.      0.2826  0.3546 -0.6653 -0.1282 -0.5791]]

最後に (3) の実装と動作確認.
(1) と (2) を繰り返して P を分解していく.
(1) で分解する回数は P の実固有値の数に一致するので予め計算しておくと便利.

# 与えられた行列の左上を I_k で水増しする関数
def pad(X, d):
    Xp = np.eye(d)
    Xp[d-np.shape(X)[0]:, d-np.shape(X)[1]:] = X
    return Xp
    
def decomp_3(P):
    d = np.shape(P)[-1]
    l = np.linalg.eigvals(P)
    dp = len(np.where(l.imag==0)[0]) # 実固有値の数 = (1) で分解する回数
    D = np.array(P)
    R = np.eye(d)
    for i in range(dp):
        Ri, Di = decomp_1(D[i:, i:])
        D[i:, i:] = Di
        R = R@pad(Ri, d)
    for i in range(dp, d, 2): # 
        Ri, Di = decomp_2(D[i:, i:])
        D[i:, i:] = Di
        R = R@pad(Ri, d)
    return R, D

R, D = decomp_3(P)
print(R)
print(D)

[[-0.1378 -0.1894  0.4886 -0.5695 -0.3763 -0.4837  0.0804]
 [-0.1928 -0.7377  0.0597  0.1701 -0.1274  0.3818  0.4735]
 [-0.7629 -0.127  -0.5124 -0.1679 -0.0257 -0.0862 -0.321 ]
 [ 0.449  -0.3475 -0.2253 -0.6084  0.1244  0.3414 -0.3532]
 [ 0.0624  0.2843 -0.1533  0.0036 -0.8647  0.3763 -0.0497]
 [-0.3952  0.3331  0.493  -0.2673  0.2526  0.5951 -0.0199]
 [ 0.0155  0.3022 -0.4217 -0.4206  0.1204 -0.0251  0.7339]]
[[-1.     -0.     -0.     -0.      0.      0.      0.    ]
 [ 0.     -0.8407 -0.5415  0.     -0.     -0.      0.    ]
 [-0.      0.5415 -0.8407 -0.      0.      0.     -0.    ]
 [-0.      0.      0.      0.2673 -0.9636 -0.     -0.    ]
 [ 0.      0.      0.      0.9636  0.2673 -0.     -0.    ]
 [ 0.      0.      0.     -0.     -0.      0.839   0.5441]
 [ 0.      0.      0.      0.      0.     -0.5441  0.839 ]]

上手に分解できました.
D の対角ブロックに先に見ていた P の固有値の実部と虚部が表れていることを確認しておこう.

この記事に素直に従うと左上ブロックは $\rho(\bar{\lambda})$ になる. $R$ の列ベクトルを適当に入れ替えれば (具体的には [実部ベクトル, 虚部ベクトル] の並びを逆にすれば) ここでの分解が導かれる. ↩
元は複素数を拡張した quaternion (四元数) という代数構造を指す. 日本語ではカタカナ語の「クオータニオン」「クォータニオン」は回転行列の代替表現という意味で使われることがほとんどなので棲み分けはできていそう. ↩
実は以上の議論は $P$ が直交行列によって対角化可能であることに頼っていて, 正規行列 ($A^tA = AA^t$ を満たす行列 $A$) であれば同様の分解が可能である. ↩
一様でなくてもいいので P = np.linalg.qr(rng.normal(size=[d, d]))[0] でもいいが. ↩