QR 分解で平均・分散共分散行列を求める

概要

OCAMI Workshop, Statistical Theories and Machine Learning Using Geometric Methods にて,

講演者「データ行列に 1 を追加した行列を QR 分解すると平均・標準偏差が出るんですが (和訳・意訳)」
ぼく「?!」

前提知識

平均と分散共分散行列が分かる程度の統計学と QR 分解が分かる程度の線型代数.
アファイン変換に関する節は付録的な位置付けなので分からなければ読み飛ばして大丈夫.

解説

手法説明

$n$ 個の $d$ 次元ベクトルデータ $x^{(j)} = [x^{(j)}_0, x^{(j)}_1, \dots, x^{(j)} _{d-1}]^t \in \mathbb{R}^d$ ($j = 0, 1, \dots, n-1$) があって $d \lt n$ の時,

$$
X := \begin{bmatrix} x^{(0)}_0 & x^{(1)}_0 & \dots & x^{(n-1)}_0 \\ x^{(0)}_1 & x^{(1)}_1 & \dots & x^{(n-1)}_1 \\ \vdots & \vdots & \vdots & \vdots \\ x^{(0)} _{d-1} & x^{(1)} _{d-1} & \dots & x^{(n-1)} _{d-1} \\ 1 & 1 & \dots & 1 \end{bmatrix} \in \mathbb{R}^{(d+1)\times n},
$$

という $(d+1)\times n$ 行列 $X$ の上三角行列による RQ 分解を考える¹.
つまり,

$$
X = RQ,
$$

であって, $R$ は $(d+1)\times(d+1)$ 上三角行列で, $Q$ は $QQ^t = I$ ($I$ は単位行列) となる $(d+1)\times n$ 行列である.

この時, $R$ は $x^{(0)}, x^{(1)}, \dots, x^{(n-1)}$ の平均ベクトル $\mu$ と分散共分散行列の上三角コレスキー因子 $U$ によって,

$$
R = \sqrt{n}\begin{bmatrix} U & \mu \\ 0 & 1 \end{bmatrix},
$$

と書ける.
つまり,

$$
\begin{gather*}
\mu := \frac{1}{n}\sum_{j=0}^{n-1}x^{(j)}, \quad \Sigma := \frac{1}{n}\sum_{j=0}^{n-1}\left(x^{(j)}-\mu\right)^2, \\
U := \mathrm{Chol}(\Sigma), \quad \text{i.e.} \ \Sigma = UU^t,
\end{gather*}
$$

であり, RQ 分解によってデータの平均・分散共分散行列を求められるということになる.

通常は RQ 分解または QR 分解の三角成分 (ここでいう $R$) について, 対角成分が全て非負であると定められることも多い.
しかしこの記事では最右下の成分が非負であることを除いて対角成分の符号は自由であっていい.

注意

この RQ 分解の計算量は $O(nd^2)$ で, 普通に平均と分散共分散行列を求める $O(nd)$ より遅い².
分散共分散の変換を目的として $U$ が欲しい場合も分散共分散行列をコレスキー分解する方が (精度はともかくとして) 速い.
逆に $1$ を追加したデータ行列の RQ 分解を求める時には有効かもしれない.

理論

$R, Q$ の $(i, j)$-成分をそれぞれ $r_{ij}, q_{ij}$ と置く.
$R$ は上三角行列なので $i \gt j$ の時 $r_{ij} = 0$ で, $Q$ の第 $i$ 行を $Q_i := [q_{i0}, q_{i1}, \dots, q_{i, n-1}]$ と置けば $QQ^t = I$ より $Q_0, Q_1, \dots, Q_d$ は互いに直交する単位ベクトルになっている³.

まず $X$ の第 $d$ 行について,

$$
\sum _{k=0}^dr _{dk}q _{kj} = r _{dd}q _{dj} = 1,
$$

より $q_{d0} = q_{d1} = \dots q_{d, n-1} = 1/r_{dd}$であり, $\|Q_d\|^2 = 1$ から $q_{dj} = 1/\sqrt{n}, r_{dd} = \sqrt{n}$ が分かる.

続いて上で求めた $Q_d = [1, 1, \dots, 1]/\sqrt{n}$ を用いると,

$$
\begin{bmatrix} \mu \\ 1 \end{bmatrix} = \frac{1}{n}X\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} = \frac{1}{\sqrt{n}}RQQ_d^t = \frac{1}{\sqrt{n}}R\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix} = \frac{1}{\sqrt{n}}\begin{bmatrix} r_{0d} \\ r_{1d} \\ \vdots \\ r_{dd} \end{bmatrix},
$$

であり, $[r_{0d}, r_{1d}, \dots, r_{d-1, d}]^t = \sqrt{n}\mu$ が得られる.

最後に, $d\times d$ 上三角行列 $R'$ を用いて,

$$
R = \begin{bmatrix} R' & \sqrt{n}\mu \\ 0 & \sqrt{n} \end{bmatrix},
$$

と書けば,

$$
XX^t = RR^t = \begin{bmatrix} R' & \sqrt{n}\mu \\ 0 & \sqrt{n} \end{bmatrix}\begin{bmatrix} R' & \sqrt{n}\mu \\ 0 & \sqrt{n} \end{bmatrix}^t = \begin{bmatrix} R'R'^t+n\mu\mu^t & n\mu \\ n\mu^t & n \end{bmatrix},
$$

となる.
この時 $i, k \leq d-1$ に対する $XX^t$ の $(i, k)$-成分は $\sum_{j=0}^{n-1}x^{(j)}_ix^{(j)}_k$ で得られるため, $R'R'^t$ の $(i, k)$-成分は $\mu = [\mu_0, \mu_1, \dots, \mu _{d-1}]^t$ という表記を用いて,

$$
\sum _{j=0}^{n-1}x^{(j)}_ix^{(j)}_k-n\mu_i\mu_k = n\left(\frac{1}{n}\sum _{j=0}^{n-1}x^{(j)}_ix^{(j)}_k-\mu_i\mu_k\right),
$$

と書かれる.
これはよく知られた共分散の公式から, $x^{(0)}, x^{(1)}, \dots, x^{(n-1)}$ の第 $i$ 成分と第 $k$ 成分に関する共分散の $n$ 倍に一致することが分かる.
よって, $R'$ は $n\Sigma$ の上三角コレスキー因子であり, $\Sigma = UU^t$ を満たす上三角行列 $U$ によって $R' = \sqrt{n}U$ と書かれる.

以上から,

$$
R = \sqrt{n}\begin{bmatrix} U & \mu \\ 0 & 1 \end{bmatrix},
$$

と書かれることが示された.

アファイン変換との関連

$d$ 次元ベクトル空間上のアファイン変換 $\varphi$ は, 線型変換 ($d$ 次正方行列) $L_\varphi$ 及び平行移動ベクトル $w_\varphi$ によって,

$$
\varphi: \mathbb{R}^d \ni x \mapsto L_\varphi x+w_\varphi \in \mathbb{R}^d,
$$

と書ける.
この時 $\varphi$ と $x \in \mathbb{R}^d$ に対して,

$$
\begin{gather*}
x^+ := \begin{bmatrix} x \\ 1 \end{bmatrix} \in \mathbb{R}^{d+1}, \\
\varphi^+: \mathbb{R}^{d+1} \ni x^+ \mapsto \begin{bmatrix} L_\varphi & w_\varphi \\ 0 & 1 \end{bmatrix} x^+ \in \mathbb{R}^{d+1}, \\
\end{gather*}
$$

という $d+1$ 次元ベクトル空間上の線型変換を考えると $\varphi^+(x^+) = (\varphi(x))^+$ が成り立ち, $d$ 次元空間上のアファイン変換は $d+1$ 次元空間上の線型変換で書き換えられる.
アファイン変換は平均と分散共分散の変換にも用いられ, 平均 $\mu$, 分散共分散行列 $\Sigma$ のデータを $\varphi$ で変換することは, 平均が $L_\varphi\mu+w_\varphi$ で分散共分散行列が $L_\varphi\Sigma L_\varphi^t$ になるように変換することに相当する.

ここで, RQ 分解 $X = RQ$ の $X$ に関する第 $j$ 列を取り出すと,

$$
\begin{bmatrix} x^{(j)}_0 \\ x^{(j)}_1 \\ \vdots \\ x^{(j)} _{d-1} \\ 1 \end{bmatrix} = \sqrt{n}\begin{bmatrix} U & \mu \\ 0 & 1 \end{bmatrix}\begin{bmatrix} q _{0j} \\ q _{1j} \\ \vdots \\ q _{d-1, j} \\ \frac{1}{\sqrt{n}} \end{bmatrix} = \begin{bmatrix} U & \mu \\ 0 & 1 \end{bmatrix}\begin{bmatrix} \sqrt{n}q _{0j} \\ \sqrt{n}q _{1j} \\ \vdots \\ \sqrt{n}q _{d-1, j} \\ 1 \end{bmatrix},
$$

であり, これは $x^{(j)}$ が $z^{(j)} := \sqrt{n}[q _{0j}, q _{1j}, \dots, q _{d-1, j}]^t$ のアファイン変換によって $x^{(j)} = Uz^{(j)}+\mu$ と書けることを意味する.
この時, 本手法の手続きで $z^{(j)}$ から生成されるデータ行列は $\sqrt{n}Q$ に一致して, その RQ 分解は $\sqrt{n}I\cdot Q$ で得られる.
従ってここまでの理論により, $z^{(j)}$ は平均 $0$ で分散共分散行列が $I$ のデータとなることが導かれ, $x^{(j)} = Uz^{(j)}+\mu$ は $z^{(j)}$ を平均 $\mu$ で分散共分散行列が $\Sigma = UU^t$ になるよう変換する式であることが分かる.
このことは $\mu, \Sigma$ が $x^{(j)}$ の平均と分散共分散行列であることに整合する.

実験

乱数生成したパラメータの正規分布に従うサンプルから直接求めた平均・分散共分散行列と RQ 分解で求めた平均・分散共分散行列の差を最大絶対値誤差で評価する.
RQ 分解は SciPy の scipy.linalg.rq によるが, $R$ の $(d, d)$-成分が負 (つまり $-\sqrt{n}$) になる可能性があるので補正を加えている.
NumPy や SciPy のコレスキー分解では下三角行列が出て来るので $U$ ではなく $\Sigma = UU^t$ を比較する必要がある.

import numpy as np
import scipy.linalg

rng = np.random.default_rng()

# サンプル数・次元
n = 1000
d = 10

# 標準正規分布に従うサンプルを変換してデータ生成
mu_true = rng.normal(size=[d, 1])
U_true = np.triu(np.exp(rng.normal(size=[d, d])))
sgm_true = U_true@U_true.T
z = rng.normal(size=[d, n])
x = mu_true+U_true@z

# サンプル平均・サンプル分散共分散行列を算出
mu_sample = np.mean(x, axis=1, keepdims=True)
sgm_sample = np.cov(x, bias=True)
# ↓だと L_sample@L_sample.T = sgm_sample なる下三角行列が出て来る
# L_sample = np.linalg.cholesky(sgm_sample)

# RQ 分解による平均・分散共分散行列を算出
# R[-1, -1] の符号で R の最終列, Q の最終行を補正する
x_aug = np.vstack([x, np.ones([1, n])])
R, Q = scipy.linalg.rq(x_aug, mode='economic')
dlt = np.sign(R[-1, -1])
R[:, -1] = dlt*R[:, -1]
Q[-1] = dlt*Q[-1]
mu_rq = R[:-1, -1:]/np.sqrt(n)
U_rq = R[:-1, :-1]/np.sqrt(n)
sgm_rq = U_rq@U_rq.T

print('max abs error of mean: {:.4e}'.format(np.max(np.abs(mu_sample-mu_rq))))
print('max abs error of covs: {:.4e}'.format(np.max(np.abs(sgm_sample-sgm_rq))))

max abs error of mean: 2.6645e-15
max abs error of covs: 4.2633e-13

参考資料

ネタ元は以下のワークショップの Mean-variance joint statistic valued in a real Siegel domain (Hiroto Inoue, Nishinippon Institute of Technology).

OCAMI Workshop - Statistical Theories and Machine Learning Using Geometric Methods

1. この設定だとタイトルと逆に RQ 分解が必要だが, タイトルでは通りの良さを考えて QR 分解とした. 転置して行と列を入れ替えれば QR 分解で良い. ↩

2. じゃあ何故この手法が提案されたかというと平均・分散共分散行列のパラメータ空間をデータ空間に作用する群と見なして事後分布に必要な積分を求める応用のためである. ↩

3. 今回はあまり影響はないが, $d \leq n-2$ も有り得るため $Q_0, Q_1, \dots, Q_d$ は $\mathbb{R}^n$ の基底になるとは限らないことに注意. このため $QQ^t$ は単位行列になるが $Q^tQ$ は常には単位行列でなく, 筆者はたまに計算を間違える. ↩