主成分分析をNIPALSアルゴリズムで解いてみる
主成分分析とは、合成変数の分散を（ある制限の元）最大化する問題を解くことである。
合成変数とは元のデータの各列を説明変数とする線形結合式で、この合成変数の分散が最大になるときの線形結合式の係数が固有ベクトルであり、そのときの合成変数を主成分と呼ぶ。
データ行列をX（n×p）とする。n > pで、ランクr = pとしておく。
X = [x1 x2 … xp] と書く。x1, x2, …, xp それぞれはXを構成する変数で行数nの列ベクトルとなる。平均0、分散1に標準化されているものとする。
合成変数をsとする。行数nの列ベクトルである。
s = f_1 x1 + f_2 x2 + … + f_p xp = X f
f_1, f_2, …, f_p が各変数に対応する係数である。f = (f_1, f_2, …, f_p)’を行数pの列ベクトルとする。
sの分散（Var(s)）を求める。
Var(s) = s’s / n-1 = (X f)’ X f / n-1 = f’ X’ X f / n-1 = f’ (X’ X / n-1) f = f’ R f
x1, x2, …, xp すべてのベクトルの平均値が0なのでsの平均値も0である。
R = X’X / n-1 は相関係数行列。
Var(s) = f’R f が最大になるときのfは、行列Rに対する固有値問題を解いたときの固有ベクトルを並べた行列（V）の第一列目の固有ベクトル。（第一列目の固有ベクトル（v1）が最大固有値（λ1）に対応する。）
V（p×r）であるが、r = pと仮定しているのでV（p×p）。
V = [v1 v2 … vp]と書く。v1, v2, …, vp それぞれはVを構成する列ベクトルで行数p。
固有ベクトルv1, v2, …, vpは、固有値λ1, λ2, …,λpに対応。
sの分散（Var(s)）の最大値がλ1である。
λ1 = v1’ R v1
sの分散が最大となる線形結合式が第一主成分（t1）となる。
第p主成分までを結合した行列をT（n×p）とする。本当はT（n×r）であるが、r = pとしているのでT（n×p）。
T = [t1 t2 … tp]と書く。t1, t2, …, tp それぞれはTを構成する主成分で行数n。
t1 = X v1
前々回、スペクトル分解を用いてR = VDV’の形で分解できると書いた。（対称行列のみ可能）
これをひっくり返してD = V’RVと書くことができる。（V’V = I）
D（p×p）は、固有値λ1, λ2, …,λpを対角要素とする対角行列。本当はD（r×r）であるが、r = pとしているのでD（p×p）。
Dからλ1のみを取り出し、Vからv1のみを取り出すと上で示したように、λ1 = v1’ R v1と書くことができる。
このように、Var(s) = f’R fを最大化するところから出発した問題が、スペクトル分解を用い、λ1 = v1’ R v1と表すことができる。（Rが対象行列なのでスペクトル分解と書いたが、特異値分解でも同じこと）
上で、t1 = X v1と書いたが、すべての主成分を作成した、と考えると、T = XVと書くことができる。
これをひっくり返してX = T V’とも書くことができる。（Vは直行行列なので、V’V = VV’ = I）
Xは、T V’という形に分解できますよ、といった見方。
この場合、X（n×p）、T（n×p）でXの列数pまで主成分を作成していることになる。主成分分析の目的は次元の縮小なので、通常pより小さい次元qまでで主成分の作成はとどめる（p > q）。
この場合、X = T* V*’ + E とデータ行列Xを書くことができる。
T*（n×q）は、主成分qまでを引っ付けた行列。V*（p×q）はq番目までの固有値（λ1, λ2, …,λq）に対する固有ベクトルの行列。
E（n×p）はXがT* V*’ で説明されない部分を示す残差行列。E = X - T* V*’。
主成分分析は、スペクトル分解や特異値分解を用いると一発で解けるのだが、主成分を一つ作成し、残差行列を求め、それに対してまた主成分を一つ作成することで計算することもできる。これがNIPALSアルゴリズム（厳密には違うかもしれないが、残差を取って繰り返し処理していくという意味でNIPALS風）による計算方法である。
Xからスタート。
X
x1 x2 x3
[1,] 0.966 0.619 -0.474
[2,] 0.207 -1.009 -0.474
[3,] -0.552 -0.358 -1.684
[4,] -0.931 -0.684 0.735
[5,] 0.587 0.944 -0.474
[6,] 1.550 1.921 1.849
[7,] -0.931 0.293 -0.474
[8,] 0.587 0.293 0.735
[9,] -1.690 -1.335 -0.474
[10,] 0.207 -0.684 0.735
Xに対し特異値分解を行い、右特異ベクトル（svd(X)$v）の第一列目を第一主成分の係数（固有ベクトル）として取り出す。
v1 <- svd(X)$v[,1]
[1] -0.6254950 -0.6104604 -0.4858951
第一主成分を作成。
t1 <- X %*% v1
[,1]
[1,] -0.75178892
[2,] 0.71679136
[3,] 1.38206543
[4,] 0.64275790
[5,] -0.71312594
[6,] -3.04063179
[7,] 0.63378525
[8,] -0.90316338
[9,] 2.10236553
[10,] -0.06905545
t1 v1’ が第一主成分まででXを近似した行列。
残差（X - t1 v1’）をE1とする。
E1 <- X - t1 %*% t(v1)
x1 x2 x3
[1,] 0.49575977 0.16006263 -0.8392906
[2,] 0.65534944 -0.57142725 -0.1257146
[3,] 0.31247506 0.48569624 -1.0124612
[4,] -0.52895813 -0.29162175 1.0473129
[5,] 0.14094327 0.50866485 -0.8205044
[6,] -0.35190008 0.06481466 0.3715719
[7,] -0.53457047 0.67990081 -0.1660469
[8,] 0.02207579 -0.25834549 0.2961573
[9,] -0.37498080 -0.05158907 0.5475291
[10,] 0.16380616 -0.72615562 0.7014463
残差E1に特異値分解を行い、右特異ベクトル（svd(E1)$v）の第一列目を第二主成分の係数（固有ベクトル）として取り出す。
v2 <- svd(E1)$v[,1]
[1] -0.2845735 -0.4013624 0.8705895
第二主成分を作成。
t2 <- E1 %*% v2
[,1]
[1,] -0.93600078
[2,] -0.06659147
[3,] -1.16530044
[4,] 1.17935314
[5,] -0.95859021
[6,] 0.39761389
[7,] -0.26532069
[8,] 0.35523946
[9,] 0.60408863
[10,] 0.85550848
[t1 t2] [v1’ v2’]’ が第二主成分まででXを近似した行列。
残差（X - [t1 t2] [v1’ v2’]’）をE2とする。計算ではE1 – t2 v2’ とする。
E2 <- E1 - t2 %*% t(v2)
x1 x2 x3
[1,] 0.22939874 -0.21561292 -0.024418072
[2,] 0.63639927 -0.59815456 -0.067740753
[3,] -0.01913857 0.01798843 0.002037182
[4,] -0.19334547 0.18172628 0.020580425
[5,] -0.13184611 0.12392276 0.014034200
[6,] -0.23874970 0.22440193 0.025413424
[7,] -0.61007371 0.57341105 0.064938561
[8,] 0.12316753 -0.11576572 -0.013110419
[9,] -0.20307318 0.19086941 0.021615880
[10,] 0.40726121 -0.38278666 -0.043350429
残差E2に特異値分解を行い、右特異ベクトル（svd(E2)$v）の第一列目を第三主成分の係数（固有ベクトル）として取り出す。
v3 <- svd(E2)$v[,1]
[1] -0.72648048 0.68282230 0.07732934
第三主成分を作成。
t3 <- E2 %*% v3
[,1]
[1,] -0.31576725
[2,] -0.87600326
[3,] 0.02634424
[4,] 0.26613994
[5,] 0.18148610
[6,] 0.32863884
[7,] 0.83976615
[8,] -0.16954004
[9,] 0.27953013
[10,] -0.56059485
[t1 t2 t3] [v1’ v2’ v3’]’ が第二主成分まででXを近似した行列。（近似と言うか、Xに一致）
残差（X - [t1 t2 t3] [v1’ v2’ v3’]’）をE3とする。計算ではE2 – t3 v3’ とする。
E3 <- E2 - t3 %*% t(v3)
x1 x2 x3
[1,] 0.000000e+00 2.775558e-17 8.326673e-17
[2,] 0.000000e+00 0.000000e+00 0.000000e+00
[3,] 4.510281e-17 1.734723e-17 2.506675e-16
[4,] 0.000000e+00 -2.775558e-17 -2.393918e-16
[5,] -5.551115e-17 -6.938894e-17 5.724587e-17
[6,] 8.326673e-17 5.551115e-17 1.040834e-17
[7,] -1.110223e-16 -1.110223e-16 1.387779e-17
[8,] 1.387779e-17 2.775558e-17 -4.683753e-17
[9,] 0.000000e+00 -2.775558e-17 -2.428613e-17
[10,] -5.551115e-17 1.110223e-16 -9.714451e-17
第三主成分まで取り出すとXをすべて再現できるので、E3はもちろん0。
私は、あいまいになるほどと思っている程度で、詳しくは理解していない。
が、PLSまでたどり着きたいので記しておく。PLSでもこれと同じ方法でスコアを作成することができる。
誤りがあれば、都度修正したい。
特異値分解は、どんな行列も非常に性質の良い3つの行列（ULV’）に分解する手法である。U, L, Vの列を元の行列のランク以下までしか使用しないと、ULV’ は元の行列の最小２乗近似となる。
もやっとしか理解できていないが、多変量解析ではとにかく使いやすいと思う。
