Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
2
Help us understand the problem. What is going on with this article?
@m1t0

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

More than 1 year has passed since last update.

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
2
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
m1t0
Major: statistics/Multivariate analysis/Psychometrics/R/ (R) C++/python3

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
2
Help us understand the problem. What is going on with this article?