« 読了：Abdi (2010) PLS回帰への招待 in 2010 | メイン | 読了： Liu & Agresti (2005) 順序カテゴリカルデータの分析方法レビュー (全体的に遠慮がちな質疑応答つき) »
2018年5月 5日 (土)
Rosipal, R., Kramer, N. (2006) Overview and Recent Advances in Partial Least Squares. in Saunders C. et al. (eds) Subspace, Latent Structure and Feature Selection. Lecture Notes in Computer Science, vol 3940. 34-51.
PLS回帰にはいろいろバージョンがあるようで、イライラするので読んでみた。先日読んだAbdiさんの解説にも引用されていた。
いわく。
PLSというのはPLS回帰だけでなく、もっと幅の広い手法である。それらの共通点は、観察データとの背後に少数の潜在変数が存在すると考える点である。ふつうは2つの変数セットについて(回帰の場合だと予測子セットと反応セット)、その間の共分散を最大化するようにして、直交する潜在変数スコアを求める。
統計学者はもともとPLS回帰を無視する傾向があり、いまでもPLS回帰を統計モデルというよりアルゴリズムとして捉えていることが多い。もっとも最近では、PCRやリッジ回帰などとあわせ、continuum回帰という統一的なアプローチとして捉えられるようになった。
PLSを分類手法としてみるとFDAと密接な関連がある。またPCAと関連付けることもできる。潜在変数をつかってSVMするという提案もある[へー]。カーネル法と併用することもできる。
標本サイズを$n$とし、$N$変数データ$\mathbf{X}$と$M$変数データ$\mathbf{Y}$があるとしよう(どちらも変数の平均は0にしてある)。PLSではこれらを
$\mathbf{X} = \mathbf{TP}^T + \mathbf{E}$
$\mathbf{Y} = \mathbf{UQ}^T + \mathbf{F}$
と分解する。$\mathbf{T}$と$\mathbf{U}$は$p$個の潜在変数のスコア行列でどちらもサイズ$n \times p$。$\mathbf{P}$と$\mathbf{Q}$は負荷行列で、サイズは$N \times p$, $M \times p$。$\mathbf{E}$と$\mathbf{F}$は残差行列。
オリジナルのNIPALSアルゴリズムでは、$\mathbf{Xw}$と$\mathbf{Yc}$の標本共分散を最大化するような、長さ1のウェイト・ベクトル$\mathbf{w}, \mathbf{c}$を求める。
手順としては、まずランダムなスコアベクトル$\mathbf{u}$を用意しておき、
- $\mathbf{w} = \mathbf{X}^T \mathbf{u}/(\mathbf{u}^T \mathbf{u})$
- $||\mathbf{w}|| \to 1$ [←これがXウェイトね]
- $\mathbf{t} = \mathbf{X w}$ [←X得点]
- $\mathbf{c} = \mathbf{Y}^T \mathbf{t} / (\mathbf{t}^T \mathbf{t})$
- $||\mathbf{c}|| \to 1$ [←Yウェイト]
- $\mathbf{u} = \mathbf{Yc}$ [←Y得点]
これを収束するまで繰り返す。なお、$\mathbf{Y}$が1変数($\mathbf{y}$)の場合には$\mathbf{u}$が$\mathbf{y}$そのものになり、繰り返しは不要になる。
実はウェイトベクトル$\mathbf{w}$は、
$\mathbf{X}^T \mathbf{Y} \mathbf{Y}^T \mathbf{X w} = \lambda \mathbf{w}$
という固有値問題の第一固有ベクトルになっている。
以上が第一因子。$\mathbf{X}, \mathbf{Y}$をデフレーションし、今度は第二因子を求める。
さて、PLSにはいろいろな形式がある。それらのちがいは、スコア$\mathbf{t}, \mathbf{u}$を求めたあとで、$\mathbf{X}, \mathbf{Y}$をどうやってデフレーションさせるかという点の違いである。4つの形式について述べる。
以下、$\mathbf{p} = \mathbf{X}^T \mathbf{t} / (\mathbf{t}^T \mathbf{t})$, $\mathbf{q} = \mathbf{Y}^T \mathbf{u} / (\mathbf{u}^T \mathbf{u})$を負荷ベクトルと呼ぶ。
- PLS Mode A。Woldのオリジナルのアプローチ。因子をひとつ求めるたびに、$\mathbf{X}$から$\mathbf{t} \mathbf{p}^T$を引き、$\mathbf{Y}$から$\mathbf{u} \mathbf{q}^T$を引く。$\mathbf{X}$と$\mathbf{Y}$は対称的に扱われている。予測というより、変数セット間の関係の記述である。正準相関分析(CCA)に近い。
- PLS1, PLS2。これは回帰手法。もっとも良く使われている手法である。$\mathbf{X},\mathbf{Y}$のどちらかが1変数の場合をPLS1,両方が 2変数以上の場合をPLS2という。
以下の2点を仮定する。(1)[$\mathbf{u}$じゃなくて]$\mathbf{t}$が$\mathbf{Y}$の良い予測子になっている。(2)スコア$\mathbf{t}$と$\mathbf{u}$のあいだに
$\mathbf{U}=\mathbf{TD} + \mathbf{H}$
という関係がある($\mathbf{D}$は$p \times p$対角行列, $\mathbf{H}$は残差行列)。
[ここから逐語訳：]予測子変数と予測対象変数のあいだのこの非対称な仮定を、デフレーションの枠組みに変換する。すなわち、予測子の空間(つまり$\mathbf{X}$)のスコアベクトル$\{\mathbf{t}_i\}^p_{i=1}$が、$\mathbf{Y}$の良い予測子である、という枠組みである。このスコアベクトルをつかって$\mathbf{Y}$をデフレートする。つまり、PLSのそれ各反復において、$\mathbf{Y}$の$\mathbf{t}$に対する回帰の要素を、$\mathbf{Y}$から引く。
$\mathbf{X} = \mathbf{X} - \mathbf{tp}^T$
$\mathbf{Y} = \mathbf{Y} - \mathbf{tt}^T \mathbf{Y} / \mathbf{t}^T \mathbf{t} = \mathbf{Y} - \mathbf{tc}^T$
ただし、ここでの$c$はNIPALSのステップ4で定義された$c$, つまり標準化する前の$c$である。
このデフレーションの枠組みは、$\{\mathbf{t}_i\}^p_{i=1}$が互いに直交になることを保証する。なお、PLS1の場合は、$\mathbf{y}$のデフレーションはテクニカルにいえば不要である。
クロス積行列$\mathbf{X}^T \mathbf{Y}$の特異値は標本共分散に対応する。そのため、一度にひとつの要素を引いていくという上記のデフレーションの枠組みは、以下の興味深い特性を持つ。$i+1$番目の反復におけるデフレートされたクロス積行列$\mathbf{X}^T \mathbf{Y}$の第一特異値は、その前の$i$番目の反復における$\mathbf{X}^T \mathbf{Y}$の第二特異値と等しいか、それよりも大きい。この結果は、上で述べた固有値の関係 $\mathbf{X}^T \mathbf{Y} \mathbf{Y}^T \mathbf{Xw} = \lambda \mathbf{w}$ にも当てはまる。なぜならこの式は$\mathbf{X}^T \mathbf{Y}$の特異値分解と対応しているからだ。特に注意したいのは、PLS1とPLS2のアルゴリズムが、上の式のすべての固有値を一度に求めるやり方とは異なるという点である。 - PLS-SB。さっき出てきた固有値問題 $\mathbf{X}^T \mathbf{Y} \mathbf{Y}^T \mathbf{Xw} = \lambda \mathbf{w}$ の固有ベクトルを一気に全部求める。このやりかただと、スコア$\mathbf{t}$は必ずしも互いに直交にならない。
- SIMPLS。de Jongが提案した方法。いちいちデフレーションせず、すべてのウェイト・ベクトル$\mathbf{w}$を$\mathbf{X}$から直接求める。スコア$\mathbf{t}$は互いに直交になる。$\mathbf{Y}$が1変量だったらPLS1と同じ。
PCA, CCAとの関係について。
PCAは$\mathbf{Xw}$の分散を最大化するような$\mathbf{w}$を求める。
CCAは$\mathbf{Xw}$と$\mathbf{Yc}$の相関を最大化するような$\mathbf{w}, \mathbf{c}$を求める。これは正準リッジ分析という概念からみることもできる。いま、
$\displaystyle max_{|r|=|s|=1} \frac{cov(\mathbf{Xr}, \mathbf{Xs})^2}{([1-\gamma_x] var(\mathbf{X r}) + \gamma_x)([1-\gamma_y] var(\mathbf{Y s}) + \gamma_x)}$
という最適化問題があるとしよう。$\gamma_x, \gamma_y$は正則化項で0以上1以下。両方を0にするとCCAとなり、両方を1にするとPLSとなる。つまり、その中間(正則化CCA)もあるわけだ。また、$Y$が1次元であればリッジ回帰もこのなかにはいる[へー！]。
なお、CCAは固有値分解で一発で求めるという点ではPLS-SBに近い。
上で述べたように、PLS1・PLS2をつかって線形回帰問題を解くことができる。
$\mathbf{X} = \mathbf{TP}^T + \mathbf{E}$
$\mathbf{Y} = \mathbf{UQ}^T + \mathbf{F}$
と
$\mathbf{U}=\mathbf{TD} + \mathbf{H}$
より
$\mathbf{Y} = \mathbf{TDQ}^T + (\mathbf{HQ}^T + \mathbf{F})$
となる。$\mathbf{C}^T = \mathbf{DQ}^T$, $\mathbf{F}^* = \mathbf{HQ}^T+\mathbf{F}$とおけば
$\mathbf{Y} = \mathbf{TC}^T + \mathbf{F}^*$
となる。これは$\mathbf{Y}$を直交予測子$\mathbf{T}$を使ったOLS回帰で分解しているのと同じである。
いま、スコアベクトル$\mathbf{t}$は正規直交ベクトルで(つまり$\mathbf{T}^T \mathbf{T} = \mathbf{I}$)、行列$\mathbf{C} = \mathbf{Y}^T \mathbf{T}$のベクトル$\mathbf{c}$は長さ1になっていないとしよう。[なぜここでこういうことを云っているのか理解できていない...]
ここで
$\mathbf{T} = \mathbf{XW}(\mathbf{P}^T \mathbf{W})^{-1}$
という関係が成り立っているので[←ここがわからない... とりあえず先に進もう]、これを代入すると
$\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{F}^*$
ただし
$\mathbf{B} = \mathbf{W}(\mathbf{P}^T \mathbf{W})^{-1} \mathbf{C}^T = \mathbf{X}^T \mathbf{U}(\mathbf{T}^T \mathbf{XX}^T \mathbf{U})^{-1} \mathbf{T}^T \mathbf{Y}$
である。[←ここもわからない...泣きたくなるね]
ここからしばらくPLS1についてだけ考える。
いま回帰式
$\mathbf{y} = \mathbf{Xb} + \mathbf{e}$
があるとしよう。$\mathbf{X}$の特異値分解を
$\mathbf{X} = \mathbf{V} \mathbf{\Sigma} \mathbf{S}^T$
として、$\mathbf{\Lambda} = \mathbf{\Sigma}^2$、その要素を$\lambda_i$としよう。また、
$\mathbf{A} \equiv \mathbf{X}^T \mathbf{X} = \mathbf{S \Lambda S}^T$
$\mathbf{z} \equiv \mathbf{X}^T \mathbf{y}$
としよう[正規方程式を簡単に書くために記号を導入しているわけね]。
さて、OLS推定量$\hat{\mathbf{b}}_{OLS}$は
$arg min_{b} || \mathbf{y} - \mathbf{Xb} ||$
の解だが、正規方程式
$\mathbf{Ab} = \mathbf{z}$
の解でもある。これを解くと
$\displaystyle \hat{\mathbf{b}}_{OLS} = \sum_i^{rank(\mathbf{A})} \frac{\mathbf{v}_i^T \mathbf{y}}{\sqrt{\lambda_i}} \mathbf{s}_i$
となる。以下では
$\displaystyle \hat{\mathbf{b}}_{i} = \frac{\mathbf{v}_i^T \mathbf{y}}{\sqrt{\lambda_i}} \mathbf{s}_i$
としよう。
線形回帰の推定量の多くは、
$\mathbf{Ab} = \mathbf{z}$
の解の近似である。たとえば主成分数$p$の主成分回帰推定量は
$\hat{\mathbf{b}}_{PCR} = \sum_i^{p}\hat{\mathbf{b}}_{i}$
だし、パラメータ$\gamma$のリッジ回帰推定量は
$\hat{\mathbf{b}}_{RR} = \sum_i^{rank(\mathbf{A})} \frac{\lambda_i}{\lambda_i + \gamma} \hat{\mathbf{b}}_{i}$
である。
ではPLS回帰はどうか... [ここからconjugate gradient法だのLanezosアルゴリズムだのKrylov空間だのという謎の呪文が乱れ飛び、到底手に負えなくなる... 数ページ飛ばします]
回帰係数$\mathbf{b}$の推定量$\hat{\mathbf{b}}$のMSEをバイアス-バリアンス分解すると、OLS推定量はバイアス0で、そのかわり($\mathbf{X}$の共分散行列の固有値が小さい時に)バリアンスが大きくなっている。そこで推定量を縮小して
$\hat{\mathbf{b}}_{shr} = \sum_i f(\lambda_i) \hat{\mathbf{b}}_i$
とすることを考える。$\lambda_i$は$\mathbf{X}$の第$i$特異値の二乗、$f(\lambda_i)$はシュリンケージ・ファクター。$i$は共分散行列のランクだけ回る。この枠組みで考えると、PCRとは「$i$が$p$以下の時$f(\lambda_i)=1$, そうでないとき0」に相当し、リッジ回帰は$f(\lambda_i) = \lambda_i / (\lambda_i + \gamma)$に相当する。
実はこの縮小推定の枠組みでPLS回帰も説明できてしまう。PLSではなんと
$\displaystyle f(\lambda_i) = 1 - \prod_{j=1}^{p} \left(1-\frac{\lambda_i}{\mu_j} \right)$
なのである。[$\mu_j$というのは... 上で読み飛ばしたくだりで説明されているんだけど、
$\mathbf{L} = \mathbf{W}^T \mathbf{A} \mathbf{W}$
の第 j 固有値らしい。では$\mathbf{W}$ってなにかというと、
$\mathbf{K} = (\mathbf{z}, \mathbf{Az}, \ldots, \mathbf{A}^{p-1}\mathbf{z})$
で張られるKrylov空間というのがあって、そこからLanczosアルゴリズムというので求めた直交基底なのだそうである... しらんがな!!!!]
このシュリンケージ・ファクターは面白い性質を持っている。$i$と$p$の組み合わせによっては$f(\lambda_i) > 1$となることさえある。$f(\lambda_i)$は$\mathbf{L}$に依存するので、つまりは$\mathbf{y}$に依存する。そのため、このシュリンケージがPLS1のMSEにどう効くのかははっきりしない。
PLSを判別・分類に使う場合は... [関心はあるんだけど、疲れちゃったのでパス]
非線形PLSとは... [カーネルPLSとか。疲れちゃったのでパス]
... というわけで、中盤でPLS回帰をいろんな線形回帰の推定量のひとつとして位置付けるあたりで頭がパンクしてしまい、肝心の後半(判別・分類や非線形PLS)は断念。また体力のあるときに。
PLS回帰の解説ってすぐに反復アルゴリズムの話になっちゃうので、そういうものかと思っていたんだけど、ちゃんとこういう統計学的な議論もあるのね。正準相関分析の正則化として捉えるとか、線形回帰の縮小推定量として捉えるとか。知りませんでした。
昔からPLS回帰とPLS-SEMとの関係がいまいちわからなくて困っているんだけど、PLS-SEMについては言及がなかった。ま、言及されても理解できないかもしれないけれどさ。
先日読んだAbdiさんの解説と見比べると、説明がちょっと違っていて混乱する... この論文で言うPLS-SBとはAbdiさんのいうBPLSのことだろうか？ Woldのオリジナルの提案は、結局$X$と$Y$を対称的に扱っていたのか(この論文)、そうでないのか(Abdiさん)？
2019/12/27: 仕事の都合で原論文を読み返し、大幅に追記。でも、正直なところさっぱりわかっていない...
論文：データ解析(2018-) - 読了: Rosipal & Kramer (2006) PLSレビュー in 2006