Posted at

多変量解析メモ-正準相関分析(Canonical Correlation Analysis, CCA)-


1. 要約

 本記事では,正準相関分析のプログラムを組むための解説記事です.自分の勉強メモとしての側面が強いので少し雑です.記事の最後にRで記述したサンプルコードを載せています.


2. 正準相関分析とは

 正準相関分析(Canonical Correlation Analysis, CCA)は,Hotelling(1936)にて開発された手法で,二つの変数群のそれぞれを混ぜ合わせ,最も相関が高くなるような変数を作成する手法である.ここで混ぜ合わされてできた変数は,正準変数(canonical variable)と呼ばれ,正準変数間の相関係数は,正準相関係数(canonical correlation)と呼ばれている.例えば,以下のように,

スクリーンショット 2019-08-29 17.45.27.png

英語,数学,体育の成績の得点を持つグループと,運動能力について測定したグループのデータに対して,正準相関分析は適用される.イマイチ私には,どう使って良いのかはピンと来ていない.赤穂(2006)のように,信号データに対して適用している例もある.

 二つの群(行列)をそれらの相関が最大になるようなベクトルに変換して,統一的な指標を得る手法と考えると,他に応用できそうな気もする.


3. 一般の定式化

以下の表記を利用する.


  • $X_{rv} = (X_1,\ldots,X_p)^T$: 確率変数ベクトル

  • $Y_{rv} = (Y_1,\ldots,Y_q)^T$: 確率変数ベクトル

  • X: 群1のデータ行列($n\times p$),行が個体,列が変数を表す.列中心化されているとする.

  • Y: 群2のデータ行列($n\times q$),行が個体,列が変数を表す.列中心化されているとする.

  • A: 群1の変数を混ぜ合わせる行列($p\times m$)

  • B: 群2の変数を混ぜ合わせる行列($q\times m$)

  • $||\cdot||_F$: フロベニウスノルム.

 正準相関分析は,正準相関係数を最大化するような係数$a,b$を求める問題として定式化される.

\begin{eqnarray}

\underset{a,b}{maximize}\ &&Corr(a^TX_{rv},b^TY_{rv})\\
&=&\frac{Cov(a^TX_{rv},b^TY_{rv})}{\sqrt{V[a^TX_{rv}]}\sqrt{V[b^TY_{rv}]}}\\
&=& \frac{a^TCov(X_{rv},Y_{rv})b}{\sqrt{a^TV[X_{rv}]a}\sqrt{b^TV[Y_{rv}]b}}\\
&=& \frac{a^TC_{xy}b}{\sqrt{a^TV_{xx}a}\sqrt{b^TV_{yy}b}}\\
subject\ to\ && a^TV[X_{rv}]a = 1,\ b^TV[Y_{rv}]b = 1
\end{eqnarray}

この問題に対するラグランジュ関数は,$\lambda_x,\lambda_y \ge 0$を用いて

$$L(\lambda_x,\lambda_y,a,b) = a^TC_{xy}b - \frac{\lambda_x}{2}(a^TV_{xx}a-1) - \frac{\lambda_y}{2}(b^TV_{yy}b-1)$$

であり,a,bのそれぞれで偏微分したものを0と置いて,停留点を求めると,

\begin{eqnarray}

\frac{\partial L}{\partial a} &=& C_{xy}b - \lambda_x V_{xx}a = {\bf 0}\tag{1}\\
\frac{\partial L}{\partial b} &=& C_{xy}^Ta - \lambda_y V_{yy}b = {\bf 0}\tag{2}
\end{eqnarray}

さらに,一つ目の式,二つ目の式にそれぞれ,$a^T,b^T$をかけることで,

\begin{eqnarray}

0&=&a^TC_{xy}b - \lambda_x a^TV_{xx}a - b^TC_{xy}^Ta + \lambda_y b^TV_{yy}b\\
&=& \lambda_y b^TV_{yy}b - \lambda_x a^TV_{xx}a
\end{eqnarray}

となる.これと,制約式$b^TV_{yy}b=a^TV_{xx}a = 1$から,

$$\lambda_x = \lambda_y$$

である.以降これらをまとめて$\lambda$とおく.仮に$V_{yy}$が正則であるとすると,(2)式から,

$$b = \frac{V_{yy}^{-1}C_{xy}^Ta}{\lambda}$$

であり,これを(1)に代入すると,

$$C_{xy}V_{yy}^{-1}C_{xy}^Ta = \lambda^2 V_{xx}a$$

が得られる.これは,$Ax = \lambda Bx$の形であり,一般化固有値問題(generalized eigen problem)と呼ばれる.これのソルバーを持っているのなら,それを使えばいいのだが,もし持っていない場合は,これを通常の固有値問題に変換することが必要となる.$V_{xx}$が正定値対称行列であることに注目すると,コレスキー分解(Cholesky decomposition)が適用できる.$V_{xx}$のコレスキー分解を

$$V_{xx} = L_xL_x^T$$

とし,$a^* = L_{x}^Ta$のように変換すると,先ほどの一般化固有値問題は,

\begin{eqnarray}

&&C_{xy}V_{yy}^{-1}C_{xy}^T(L_{x}^T)^{-1}a^* = \lambda^2 L_{x}a^*\\
&\Leftrightarrow& L_{x}^{-1}C_{xy}V_{yy}^{-1}C_{xy}^T(L_{x}^T)^{-1}a^* = \lambda^2 a^*
\end{eqnarray}

のように通常の固有値問題に帰着される.この問題の解を$a^* = a_{ans}$とすると,正準相関係数$a$の解は,

$$\widehat{a} = (L_x^T)^{-1}a_{ans}$$

で表される.$b$も同様のステップで求められる.


4. 別の定式化(SVDを用いる)

正準相関分析は,以下の最小二乗法基準を最小化することでもパラメータを推定できる.

\begin{eqnarray}

\underset{A,B}{minimize}&&\|XA-YB\|_F^2 = tr(XA-YB)^T(XA-YB)\\
subject\ to\ &&\frac{1}{n}A^TX^TXA =I_m,\ \frac{1}{n}B^TY^TYB=I_m
\end{eqnarray}

この定式化のざっくりとしたイメージは,合成して出来た行列が似ている(=ダイバージェンスが小さい)と,正準相関係数も高いというものである.損失関数は,

\begin{eqnarray}

\|XA-YB\|_F^2 &=& tr(XA-YB)^T(XA-YB) \\
&=& tr(A^TX^TXA - A^TX^TYB + B^TY^TYB) \\
&=& tr(nI_m - A^TX^TYB + nI_m)\\
&=& 2nm - trA^TX^TYB
\end{eqnarray}

のように展開されることから,損失関数の最小化は,$trA^TX^TYB$の制約付き最大化問題に帰着する.ここで,

\begin{eqnarray}

A^TX^TYB &\propto& A^T(n^{-1}X^TY)B\\
&=& A^TC_{xy}B\\
&=& A^TV_{xx}^{1/2}V_{xx}^{-1/2}C_{xy}V_{yy}^{-1/2}V_{yy}^{1/2}B\\
&=& A^TV_{xx}^{1/2}R_{xy}V_{yy}^{1/2}B\\
&&(R_{xy} := V_{xx}^{-1/2}C_{xy}V_{yy}^{-1/2})\\
&=& A_*^TR_{xy}B_*\\
&&(A_*: = V_{xx}^{1/2}A,\ B_*:=V_{yy}^{1/2}B)
\end{eqnarray}

であることから,一般的な正準相関分析の定式化に対応していることがわかる.

 制約を満たして,$trA_* ^TR_{xy}B_*$を最大化するA,Bは,特異値分解$$R_{xy} = KDL^T$$を用いて,

\begin{eqnarray}

A &=& V_{xx}^{-1/2}K_t\\
B &=& V_{yy}^{-1/2}L_t
\end{eqnarray}

と表される.ここでtは分析者があらかじめ設定する正準変数の個数である.また,特異値分解の結果の対角行列$D$の対角要素が正準相関係数である.


5. サンプルコード(R言語)

 特異値分解を用いる定式化のプログラムを書きました.データと一部コードはhttps://stats.idre.ucla.edu/r/dae/canonical-correlation-analysis/を参考にさせていただいています.参考にした部分は"="が"<-"になっています.


R,cca.R

CCA.my = function(X,Y,ndim){

#suppose that nx = ny
nx = nrow(X); ny = nrow(Y)
px = ncol(X); py = ncol(Y)
for(i in 1:px){
X[,i] = X[,i] - mean(X[,i])
}
for(i in 1:py){
Y[,i] = Y[,i] - mean(Y[,i])
}
Cxx = (1/nx) * t(X) %*% X
Cyy = (1/ny) * t(Y) %*% Y
Cxy = (1/nx) * t(X) %*% Y
evd.Cxx = eigen(Cxx)
Cxx.sqrt.inv = evd.Cxx$vectors %*% diag(1/sqrt(evd.Cxx$values)) %*% t(evd.Cxx$vectors)
evd.Cyy = eigen(Cyy)
Cyy.sqrt.inv = evd.Cyy$vectors %*% diag(1/sqrt(evd.Cyy$values)) %*% t(evd.Cyy$vectors)
R = Cxx.sqrt.inv %*% Cxy %*% Cyy.sqrt.inv
svd = svd(R)
A = Cxx.sqrt.inv %*% svd$u[,1:ndim]
B = Cyy.sqrt.inv %*% svd$v[,1:ndim]
return(list(canonical.corr = svd$d[1:ndim], canonical.coef.A = A, canonical.coef.B = B))
}

############################################################################

mm <- read.csv("https://stats.idre.ucla.edu/stat/data/mmreg.csv")
colnames(mm) <- c("Control", "Concept", "Motivation", "Read", "Write", "Math",
"Science", "Sex")
psych <- mm[, 1:3]
acad <- mm[, 4:8]

library(CCA)#既存のパッケージ
cc1 <- cc(psych, acad)
cc1$cor
cc1$xcoef
cc1$ycoef

CCA.my(as.matrix(psych), as.matrix(acad),3)



参考文献


  • Hoteling(1936).Relation between two sets of variates, Biometrika, 28, 321-277

  • Hardon, Szedmak, and Shawm-Taylor(2004). Canonical Correlation Analysis: An Overview with Application to Learning Methods, Neural Computation 16, 2639-2664

  • 赤穂(2006). 正準相関分析入門 —複数種類の観測からの共通情報抽出法—, 日本神経回路学会誌, 20(2), 62-72

  • 山下(2012). 正準相関分析の相関平方和最大化による定式化と構造行列の回転, 行動計量学, 39, 1-9