数学から見た対応分析

最近対応分析（コレスポンデンス分析）を学んだのですが，数学的に綺麗でない説明ばかりで分かりにくかったので，私なりに紙とペンで計算してみました．
"数学的に綺麗" というゴールにはまだ遠いですが，行列を用いてすっきり計算できたのではないでしょうか．

はじめに

対応分析とは，$m \times n$ の分割表が与えられたとき，その行と列がともに $p$ 個の "潜在変数" で説明されると考えて，$p$ 次元空間にプロットする手法を言います．
ある行と列が近くにプロットされている場合，これらの間に相関があることになります．

対応分析の統計的側面（どのような検定を行うのか等）や実例については，分かりやすい解説がたくさん溢れているのでそちらにお任せするとして，ここでは数学的な under the hood を調べていこうと思います．

行や列を表す大きな空間 $\mathbb{R}^m, \mathbb{R}^n$ をより小さな次元の空間 $\mathbb{R}^p$ へ変換するわけですが，その際大事になるのはこの変換が「何を保つ」変換なのかです．
対応分析では行の間の $\chi^2$ 距離と列の間の $\chi^2$ 距離を保とうとします．

もう少し一般に "$\chi^2$ 内積" を考え，これを保つ変換，という制約を課します．
すると非常に自然な流れで，対応分析において重要な諸量が導かれます．

※"内積を保つ" は不正確な表現である可能性が高いので，適切な解釈を考え中です．(2022-08-20)

記法

数ベクトル空間 $\mathbb{R}^n$ の標準基底を，列ベクトルと思うときは $\{ e_1, \dotsc, e_n \}$ として，行ベクトルと思うときは $\{ e^1, \dotsc, e^n \}$ と書きます．
また，列ベクトルの成分を上付き添字で，行ベクトルの成分を下付き添字で表します．たとえば $v = \sum v^i e_i$，$w = \sum w_i e^i$ などです．

行列の成分は $A = (a^i_j)$ のように上下の添字で表します．このとき $a^i := e^i A = (a^i_1, \dotsc, a^i_n)$ は $A$ の $i$ 番目の行ベクトルで，$a_j := A e_j = (a^1_j, \dotsc, a^m_j)^{\mathrm T}$ は $j$ 番目の列ベクトルになります．また，$a^i_j = e^i A e_j$ が成り立ちます．

以下，原則として添字 $i$ は上付き（"行" を表す），$j$ は下付き（"列" を表す）とします．

入力

頻度表として $m \times n$ 行列 $A = (a^i_j) \in \mathrm{Mat} (m \times n, \mathbb{R}_{\geqslant 0})$ が与えられます．
ただし，"頻度 $0$" のカテゴリーは除いておきます．つまり，各列ベクトル $A e_j$ および行ベクトル $e^i A$ は $\neq 0$ とします．

平均と重み

以下では "行カテゴリー" $a^i := e^i A = (a^i_1, \dotsc, a^i_n) \in \mathbb{R}^n$ を中心に考えますが，"列カテゴリー" $a_j := A e_j \in \mathbb{R}^m$ を中心として考えても構いません．

"全頻度" を $A_0 := \sum_{i, j} a^i_j$ として，この全体に対する相対頻度行列を $P := (1 / A_0) A$ とします．

さらに，全体のうちカテゴリー $i$ の占める割合を $\mu^i := \sum_j a^i_j / A_0$ として，カテゴリー $i$ におけるカテゴリー $j$ の割合を $p^i_j := a^i_j / \mu^i A_0$ とします．
すると，行ベクトル $p^i := (p^i_1, \dotsc, p^i_n) = (1 / \mu^i A_0) a^i \in \mathbb{R}^n$ を得ます．
この $\mu^i$ をカテゴリー $i$ の "重み" あるいは質量 (mass) と呼び，各 $p^i$ をプロファイル (profile) と呼びます．

また，カテゴリー $j$ の周辺分布を $\mu_j := \sum_i p^i_j \mu^i = \sum_i a^i_j / A_0$ と置くと，平均プロファイル $\mu := (\mu_1, \dotsc, \mu_n) \in \mathbb{R}^n$ を得ます．

便宜上，重みのなす行列を２つ定義しておきます：

$$\begin{align*} M_{\mathrm R} &:= \mathrm{diag} (\mu_1, \dotsc, \mu_n) \in \mathrm{Mat} (n, \mathbb{R}_{\gt 0}), \\ M_{\mathrm C} &:= \mathrm{diag} (\mu^1, \dotsc, \mu^m) \in \mathrm{Mat} (m, \mathbb{R}_{\gt 0}) . \end{align*}$$

$\chi^2$ 距離と慣性

プロファイル間の距離を，$\chi^2$ 距離 $d_\chi (p^i, p^k)^2 := \sum_j (p^i_j - p^k_j)^2 / \mu_j$ で測ることにします．
この距離は，内積 $\langle p^i \mid p^k \rangle := p^i M_{\mathrm R}^{-1} p^{k \mathrm T}$ から定まる距離と一致します．

対応分析では，この $\chi^2$ 距離を用いて定義される慣性モーメント $I := \sum_i \mu^i d_\chi (p^i, \mu)^2$ を不変にしながら，より少ないパラメータ空間 $\mathbb{R}^p$ へ変換することを考えます．
ここで，プロファイル $p^i$ と平均プロファイル $\mu$ との距離は

$$d_\chi (p^i, \mu)^2 = \langle p^i - \mu \mid p^i - \mu \rangle = (p^i - \mu) M_{\mathrm R}^{-1} (p^i - \mu)^{\mathrm T}$$

で与えられるため，より一般に，次で定まる $m$ 次行列 $X$ を使います：

$$X := \Delta M_{\mathrm R}^{-1} \Delta^{\mathrm T}, \quad \Delta := \begin{pmatrix} p^1 - \mu \\ \vdots \\ p^m - \mu \end{pmatrix} \in \mathrm{Mat} (m \times n, \mathbb{R}).$$

この行列の $(i, j)$ 成分は，$e^i X e_j = \langle p^i - \mu \mid p^j - \mu \rangle$ となります．
したがって，慣性モーメントが

$$\begin{align*} I &= \sum_i \mu^i e^i X e_i \\ &= \sum_i e^i M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2} e_i \\ &= \mathrm{tr} \left( M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2} \right) \end{align*}$$

で計算できます．

ついでに慣性中心のノルムも計算しておきましょう．慣性中心 $v = \sum \mu^i e^i$ について

$$\begin{align*} v \Delta = \sum_i \mu^i e^i \Delta = \sum_i \mu^i (p^i - \mu) \end{align*}$$

が成り立ちますが，定義より $\sum_i \mu^i p^i = \sum_i a^i / A_0 = \mu$ かつ $\sum_i \mu^i = 1$ となるので，$v \Delta = 0$ が導かれます．
とくに慣性中心のノルムは $v X v^{\mathrm T} = (v \Delta) M_{\mathrm R}^{-1} (\Delta^{\mathrm T} v^{\mathrm T}) = 0$ になります．

標準化残差

上の計算を見て見ると，慣性モーメントよりも少し一般的に，$\chi^2$ 内積を与える行列 $X$（または $M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2}$）を考えればよいことが分かります．

さて，定義より $e^i X e_i = d_\chi (p^i, \mu)^2 \geqslant 0$ なので，$M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2}$ は半正定値行列であり，したがって行列 $S \in \mathrm{Mat} (m \times n, \mathrm{R})$ による分解 $M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2} = S S^{\mathrm T}$ を持ちます．
特に，慣性モーメントは $I = \mathrm{tr} (S S^{\mathrm T}) = \sum_{i, j} (s^i_j)^2$ で求まります．

標準化残差 (standardized residuals) $S$ を具体的に書いてみましょう．$X = \Delta M_{\mathrm R}^{-1} \Delta^{\mathrm T}$ を "真っ二つ" に分けて $S = M_{\mathrm C}^{1/2} \Delta M_{\mathrm R}^{-1/2}$ と置けば良いのですが，この表式は a priori には "行か列か" の選択に依存しているので，少し気持ち悪いです．
そこで

$$W := \mu \begin{pmatrix} \mu^1 \\ \vdots \\ \mu^m \end{pmatrix} = \begin{pmatrix} \mu^1 \mu_1 & \cdots & \mu^1 \mu_n \\ \vdots & & \vdots \\ \mu^m \mu_1 & \cdots & \mu^m \mu_n \end{pmatrix} \in \mathrm{Mat} (m \times n, \mathbb{R})$$

という行列を導入して，$P - W$ を考えてみます．
$P - W$ の第 $(i, j)$ 成分が $(a^i_j / A_0) - \mu^i \mu_j = \mu^i (p^i_j - \mu_j)$ なので，$P - W = M_{\mathrm C} \Delta$ が成り立っています．よって $S = M_{\mathrm C}^{-1/2} (P - W) M_{\mathrm R}^{-1/2}$ が成り立ち，"行と列" の扱いに関して対称な表式を得ます．

とくに列プロファイルの間の $\chi^2$ 距離 $d_\chi^{\mathrm C} (p_j, p_\ell)^2$ について，上と同じように "$\chi^2$ 内積" を表す行列 $X_{\mathrm C} = \Delta_{\mathrm C} M_{\mathrm C}^{-1} \Delta_{\mathrm C}^{\mathrm T}$ を定めたときに，$M_{\mathrm R}^{1/2} X_{\mathrm C} M_{\mathrm R}^{1/2} = S^{\mathrm T} S$ が成り立ちます．

特異値分解

※この節の内容は不確かな部分があります．(2022-08-20)

証明は省きますが，実行列 $S$ に対して特異値分解 (singular value decomposition, SVD) $S = U \Sigma V^{\mathrm T}$ が可能です．
ここで $U \in O (m \times p)$ と $V \in O (n \times p)$ は半直交行列 (semi-orthogonal matrices)

$$U^{\mathrm T} U = 1, \quad V^{\mathrm T} V = 1$$

であり，$p = \mathrm{rank}\, S$ は $S$ の階数，$\Sigma = \mathrm{diag} ( \sigma_1, \dotsc, \sigma_p ) \in \mathrm{Mat} (p, \mathbb{R}_{\gt 0})$ は対角行列です．

※元々 $U U^{\mathrm T} = 1$ などど書いていましたが，そのような等式はあり得ません．（両辺の rank を比較すれば良い．）(2022-08-20)

今 $p \leqslant \min \{ m, n \} - 1$ に注意しましょう．実際，行プロファイルの慣性中心 $v = \sum_i \mu^i e^i$ に対して $v \Delta = 0$ が成り立っていたので，$\mathrm{rank}\, S = \mathrm{rank}\, \Delta \leqslant m - 1$ が分かります．
同様に $\mathrm{rank}\, S \leqslant n - 1$ なので，結局 $p \leqslant \min \{ m, n \} - 1$ が従いました．

特異値分解の意味を解釈するために，行カテゴリーの $p$ 次元表示 $F := M_{\mathrm C}^{-1/2} S V$ を考えます．各行ベクトル $f^i = (f^i_1, \dotsc, f^i_p) := e^i F$ をプロファイル $p^i$ の $p$ 次元空間 $\mathbb{R}^p$ での座標と思うと，Euclid 内積（$L^2$ 内積）$f^i f^{j \mathrm T}$ が $\chi^2$ 内積 $\langle p^i - \mu \mid p^j - \mu \rangle$ と一致します：

$$\begin{align*} f^i f^{j \mathrm T} &= e^i F F^T e_j \\ &= e^i M_{\mathrm C}^{-1/2} S V V^{\mathrm T} S^{\mathrm T} M_{\mathrm C}^{-1/2} e_j \\ &= e^i X e_j = \langle p^i - \mu \mid p^j - \mu \rangle . \end{align*}$$

とくに慣性中心 $v = \sum p^i e^i$ の座標表示 $w = \sum p^i f^i$ について $w w^{\mathrm T} = v X v^{\mathrm T} = 0$ が成り立つので，$w = 0$，つまり慣性中心が原点に来ることになります．

同様に列カテゴリーの $p$ 次元表示 $F_{\mathrm C} := U^{\mathrm T} S M_{\mathrm R}^{-1/2}$ を取ると，行プロファイルの $p$ 次元座標 $f_j^{\mathrm C} := F_{\mathrm C} e_j$ の $L^2$ 内積が $\chi^2$ 内積と一致します：$f_i^{\mathrm C\, \mathrm T} f_j^{\mathrm C} = e^i X_{\mathrm C} e_j$．

以上の計算から，半直交行列 $U, V$ は "$\mathbb{R}^m$ や $\mathbb{R}^n$ から $\mathbb{R}^p$ への写像で，$\mathbb{R}^m, \mathbb{R}^n$ における内積（$X, X_{\mathrm C}$）を保つ写像" であると言うことができます．

固有値と固有ベクトル

特異値分解において，$\sigma_1, \dotsc, \sigma_p$ のことを $S$ の特異値 (singular values)，特異値を $S v = \sigma_i u$ のように与えるベクトル $0 \neq u \in \mathbb{R}^m$ と $0 \neq v \in \mathbb{R}^n$ を特異ベクトル (singular vectors) と呼びます．
これら（or 2乗）は $S^{\mathrm T} S$ や $S S^{\mathrm T}$ の固有値・固有ベクトルとなるため，対応分析では単に固有値・固有ベクトルと呼ぶことも多いです．

特異値の2乗の和 $\sum_i \sigma_i^2$ が $\mathrm{tr} \, \Sigma^2$ であることを踏まえると，慣性モーメントが

$$\begin{align*} I &= \mathrm{tr} (M_{\mathrm C}^{1/2} X M_{\mathrm C}^{1/2}) \\ &= \mathrm{tr} \, S S^{\mathrm T} \\ &= \mathrm{tr} \, \Sigma \Sigma = \sum_i \sigma_i^2 \end{align*}$$

と計算できます．

また，各特異値は $\sigma_k^2 = \sum_i \mu^i (f^i_k)^2$ を満たします．
それを確認するために，まず $F = M_{\mathrm C}^{-1/2} S V = M_{\mathrm C}^{-1/2} U \Sigma$ と変形し，$F$ の列ベクトル $f_j := F e_j$ を考えます．すると

$$\begin{align*} \sum_i \mu^i (f^i_k)^2 &= f_k^{\mathrm T} M_{\mathrm C} f_k \\ &= e^k F^{\mathrm T} M_{\mathrm C} F e_k \\ &= e^k \Sigma^2 e_k = \sigma_k^2 \end{align*}$$

が分かります．

おわりに

以上，対応分析の内部機構を計算してきました．成分表示による計算とは違い，行列だとかなり見通しのよい計算になりますね．
"$\chi^2$ 内積を保つように $p$ 次元で表現する" という目標にしたがって，すべてのステップが自然に流れていきました．

次なるゴールは "座標表示（基底の取り方）に依らない線形不変な" 定式化ですが，ここまででも十分でしょう．