部分的最小二乗回帰の話をしよう~導出とアルゴリズム~

この記事でする事

部分的最小二乗回帰の事をここではPLSと略します。

• 単純な逆問題の解法と問題点
• PLSとは何か
• PLSの使い道
• PLSの導出
• PLSの実装
• PLSの過学習と対策(まだ書いていないです)

これはこの記事を書くにあたって書いたコードです。

部分的最小二乗回帰の話をしよう。

部分的最小二乗回帰は逆問題の欠点の一部を解決した解法の一つです。

これには幾通りもの導出方法とアルゴリズムがあるから何から話したものか。多分、一番有名なアルゴリズムはnipals、次点で特異値分解。

この記事は数学については詰めていませんし、逆行列みたいな基礎を書いていたりします。何故かって？神は言っている。「全ての数学弱者をも救え」と。

単純な逆問題の解法と問題点

ここ、分る人は「問題点」までとばして下さい。

逆問題って何？

問1. 鶴と亀が居ます。個体の合計は8匹、足の合計は22本です。鶴と亀、それぞれ何匹居るかを答えよ。

解1. 鶴5匹 亀3匹

この問1を鶴亀算🐢といいます。いにしえの受験戦士の一部が慣れ親しんだ演算ですが奥が深いです。下記が似たような感じであつかえます。

• 鶴亀算
• 連立方程式の解を求める事
• 行列の割り算をする事

では、鶴亀算を連立方程式にして割り算にしましょう。

$$A= \left(\begin{array}{rr} 2本 & 4本 \\
1体 & 1体 \end{array}\right) $$

$$X= \left(\begin{array}{r} x匹 \\
y匹 \end{array}\right) $$

$$Y=\left(\begin{array}{r} a本 \\
b体 \end{array}\right)$$

$$AX=Y$$
逆行列や疑似逆行列を用いて
$$X = A^{\dagger}Y$$

はい、割り算ですね。

割り算の問題点

問1を再掲します。

問1. 鶴と亀が居ます。個体の合計は8匹、足の合計は22本です。鶴と亀、それぞれ何匹居るかを答えよ。

では問2を出します。

問2. 「鶴と亀と猫が居ます。個体の合計は8匹、足の合計は22本です。鶴と亀と猫、それぞれ何匹居るかを答えよ。」…こんな問題で大丈夫か？大丈夫でないならその理由を述べよ。

解2.大丈夫じゃない、問題だ！解が無限に存在する不良設定問題だからだ！

具体的には亀と猫の比率の自由度が高くなりすぎて解を導けません。例えば「猫1匹、亀2匹」が答えだったとしても「猫103匹、亀-100匹」という解を弾き出せてしまうのです。つまり、根本的にはあまりに似た者どうしが方程式の中に入ると不都合なのです。
ちなみに、疑似逆行列を使うとそれぞれ5匹、1.5匹、1.5匹になります。

そんなの当たり前じゃねーかと思われそうですが、もっとファジーな近しいものがある時に困ります。
別の例えで言うと…脳と成績に関するものにしてみましょうか。僕が高校生の頃に履修した科目は以下の通りです。

「英語, 現代文, 古文, 漢文, 算数, 物理, 化学, 生物, 地理, 日本史, 保健体育, 体育」
さて、これらの成績と頭部MRIの結果を比べるとしましょう。頭部MRIは三次元空間のなかの座標によって、水分子の磁気の変化を記録する装置です。つまり、頭部MRIの出力結果は下記のようになりますね。

• X軸
• Y軸
• Z軸
• 磁場

これでMRIの出力が出たとしましょう。「凄いね、算数している時は脳のここが動くんだね！」とかなって面白そう。でも算数と物理は似ていますよね？眺める分にはいいですが、方程式を解くとなると算数と物理の比率がおかしな結果になりかねません。この時、算数と物理が亀と猫の関係に近くなっています。

むりやり割り算する

大抵の場合は先の鶴亀算のようにうまくはいきません。鶴亀猫算のような計算では綺麗に計算できないのです。鶴亀算は特殊な場合なのです。無限に答えがでてしまう場合はもっとも確からしい答えを探します。

下記の条件を満す時に不良設定問題と言います。

• 解が存在する
• 解が一意である
• パラメタが連続的に変化したとき、解も連続的に変化する

鶴亀カブト虫算には解が無限に存在するため、不良設定問題と言えます。逆に解が無い場合もあります。その場合は一番良い解っぽいものを求める事になります。
https://ja.wikipedia.org/wiki/%E8%89%AF%E8%A8%AD%E5%AE%9A%E5%95%8F%E9%A1%8C

$Y = AX$

ここで、Xの2乗は小さい方がいいぞ、とするのが最小ノルム解です。そして、

$A^\dagger Y = X$

となる$A^\dagger$が疑似逆行列えす。でも疑似逆行列と言われてもアレなので導出します。

ラグランジュの未定乗数法とは、極大、極小を求める時に使えるテクニックですね。先に微分して云々とか言ってたのと似た手法です。

https://manabitimes.jp/math/879

ここの$\lambda$はラグランジュの未定乗数法でよく出てくる適当な値です。

$L=X^TX - \lambda (Y-AX)$
$\frac{dL}{dX} = X + A^T\lambda^T=0$
$X = -A^T\lambda^T$
$Y=AX = -AA^T\lambda^T$
$\lambda^T = -(AA^T)^{-1}Y$
$X = A^T(AA^T)^{-1}Y$

というわけで、この場合は$A^T(AA^T)^{-1}$が割り算ですね。疑似逆行列の一部です。

今回は解が無限に存在するケースなので最もたしからしい解を探しますが、解が存在できないケースもあります。鶴亀猫は無限の解、鶴亀算は一意の解、そして鶴算は解なしです。解なしの場合も最も確からしい解を求めます。

ってなりますけど、根本的に「四足の子」とか「理系能力」みたいなのが複数あるのが悪いんですよね。

潜在変数

問3. 亀と猫みたいに似た者が居ても、大丈夫で問題なくなるようにする方法を述べよ

解3. 「同じようなもの」は一つにまとめる。

まとめるために使う変数を潜在変数といいます。具体的に、方程式をこんな風に分解してみましょう。

$$\left(\begin{array}{r}
科目ごとの成績
\end{array}\right)
\leftarrow
\left(\begin{array}{r}理解力\\
読解力\\
暗記力\\
忍耐力\\
体力
\end{array}\right)
\Leftrightarrow
\left(\begin {array}{r}
前頭葉\\
海馬\\
側頭葉\\
扁桃体
\end{array}\right)
\rightarrow
\left(\begin {array}{r}
X軸\\
Y軸\\
Z軸\\
磁場
\end{array}\right) $$

どうですか？この真ん中の二つのパラメタ、面白そうじゃあないですか？これらがここで言う潜在変数です。

よっしゃ!じゃあ潜在変数を求めるぞ!ってなりますよね?ただこれ、行列の割り算だけじゃちょっと厳しいんですよねー。この割り算の中には何も潜在変数なんて出てきません。

$$X = A^{\dagger}Y$$

潜在変数なんて無いじゃないかと思いましたか？無いなら作ればいいんです。大丈夫だ、問題ない。潜在変数はありまぁす。

前提

ここで、潜在変数を$Z$とします。X関連の$Z$を$Z_x$、Y関連の$Z$を$Z_y$とします。これらは用語では「スコア」と呼ばれるそうです。また、前提として、ここで$X$と$Y$はあらかじめスケールしておきます。具体的には平均値を引いた上で標準偏差で割り算します。
で、$w_x$と$w_y$をそれぞれ潜在変数に変換するための重み行列としておきます。

• $X w_x = Z_x$
• $Y w_y = Z_y$

また、下記のように$X_L$と$Y_L$を作ります。

• $X = Z_x X_L + E$
• $Y = Z_y Y_L+ F$

ただし、$E$と$F$は誤差項です。本来はこの式が先だったり、$Z_y$は$U$と表記されていたりしますが、$Z_y$みたいな表記の方がここではわかりやすいかなって。$X_L$と$Y_L$は後で使います。

さて、こうなる$w_x$, $w_y$を求めたいですね。でもその前に$w_x$と$w_y$の前提を確認します。$w_x^T w_x=I$で$w_y^Tw_y=I$とします。あくまで「重み」とか「比率」的ななにかですからね。どんな風に求めましょうか？部分的最小二乗回帰では、$Z_x$と$Z_y$の共分散が大きくなる解を「一番良い解」とします。先程の数式を再掲します。

$$\left(\begin{array}{r}
科目ごとの成績
\end{array}\right)
\leftarrow
\left(\begin{array}{r}理解力\\
読解力\\
暗記力\\
忍耐力\\
体力
\end{array}\right)
\Leftrightarrow
\left(\begin {array}{r}
前頭葉\\
海馬\\
側頭葉\\
扁桃体
\end{array}\right)
\rightarrow
\left(\begin {array}{r}
X軸\\
Y軸\\
Z軸\\
磁場
\end{array}\right) $$

この例えで言うと真ん中の太矢印の間の関係性が強いものこそが「一番良い解」であるわけです。
先にも言った通り、これには複数のアルゴリズムがありますが、本記事では二通り紹介します。

問4. 潜在変数$Z_x$と$Z_y$の共分散が最大になるように潜在変数$Z_x$と$Z_y$を求めよ。ただし、重み行列のノルムは1とする。

nipals

問4.1 部分的最小二乗回帰におけるnipalsアルゴリズムを述べ、既存のソフトウェアと同様に動くように実装せよ。なお、求める潜在変数は一部分でよいものとする。

ここではnipalsという方法を見ます。後で特異値分解の方法を見ます。
潜在変数$Z_x$と$Z_y$の関連を最大にする式をまず作ります。

$$Cov(Z_x, Z_y) = Z_y^T Z_x = Z_y^T X w_x$$

ここで、$w_x^T w_x = w_y^T w_y = 1$、にします。「え？」と思われるかもしれませんが、ここで重み行列は1行のものにしておきます。一行ずつ$w_x$を求めて、最終的に全部求めて行列にしたらよいのです。

求め方

さて、ここで一つ困った事があります。決めなきゃいけない変数が多すぎるのです。通常の高校数学では方程式の解は一つずつしか解けないのです。

そこで、$Z_y$を適当な値にします。

$$Z_y = 今俺が適当に決めた数$$

やばい！こんな仮定で大丈夫か？
大丈夫だ、問題ない。

この仮定のもと、一番良い$w_x$を求めます。

では$Cov(Z_x, Z_y)$に対してラグランジュの未定乗数法を使います。

$$\frac{dL}{dw_x} = Z_y^T X - 2\lambda w_x = 0$$ $$\frac{Z_y^T X}{2\lambda} = w_x$$ ここで、$w_x$はノルム(敢て数学者に殺される表現をすると「大きさ」の事)が1の傾きだけの存在です。また、$w_y$が1行なので$Z_y$も1行ですね。つまり、$2\lambda$は適当に$w_x$のノルムを1にするような数となります。

$$w_x' = {Z_y^T X}$$
$$w_x = \frac{w_x'}{||w_x'||}$$

この$w_x$は「$Z_y=今俺が適当に決めた数$」というヤバい条件下での最良の$w_x$です。

前提がヤバくても、その中で最良を求める事はできます。なので$X w_x = Z_x$より最良の$Z_x$も計算できます。では、ここで再びラグランジュの未定乗数法で今度はXとYを逆にして、最良の$w_y$を求めてみましょう。
$w_y' = {Z_x^T X}$
さきほど、一番良い$w_x$から$Z_x$を計算できていますから、これにより一番良い$Z_y$が分ります。では最後に$w_y$から$Z_y$を計算しましょう。

$$Z_y = X w_x \ne 俺がさっき適当に決めたZ_y$$

解は違っている、ここで計算を終える定めではないぞ！

巻き戻す

(一つ前のはじまりまで戻します)
そこで、$Z_y$を適当な値にします。

$$Z_y = 今俺が適当に決めた数$$
やばい！こんな仮定で大丈夫か？

そんな訳ない。一番良い数を頼む！

$$Z_y = 今俺が現時点で知っている一番いい数 = 先程計算したZ_y$$

この仮定のもと、一番良い$w_x$を求めます。
(以下、これを沢山やる)
こうして「俺が適当に決めた数」から「その中で一番良い結果」を求めていき、少しずつ正解に近くなっていくのがnipalsです。

pythonで書くとおおよそこんなコードになります。このコードは試しに自分で書いてみたコードです。インデントが変ですが、コードの途中を抜いてきたものですから勘弁してください。

    def _fit1(self, x, y):
        y_score = self.y[:, 0:1]
        for n in range(self.iter_num):
            # x_weight = ((y_score.T @ self.x) / (y_score.T @ y_score)).T
            x_weight = (y_score.T @ self.x).T
            x_weight /= np.linalg.norm(x_weight)
            x_score = x @ x_weight
            # y_weight = (self.y.T @ x_score) / (x_score.T @ x_score)
            y_weight = self.y.T @ x_score
            y_weight /= np.linalg.norm(y_weight)
            y_score = y @ y_weight
        self._flip_weight([x_weight, y_weight, x_score, y_score])
        x_load = (x_score.T @ self.x) / (x_score.T @ x_score)
        y_load = (y_score.T @ self.y) / (y_score.T @ y_score)

さて、一つだけ問題があります。$Z_y$は1行の行列でした。そんな一部だけ取ってきた行列を解として大丈夫か？実は全然大丈夫じゃないので後述します。

特異値分解を使う方法

問4.2 特異値分解を用いて部分的最小二乗回帰のアルゴリズムを述べ、既存のソフトウェアと同様に動くように実装せよ。なお、求める潜在変数は一部分でよいものとする。

ここでは特異値分解を使う方法を見ます。

$$Cov(Z_x, Z_y)=Cov(X w_x, Y w_y)=w_x^T X^T Y w_y$$

ここで、$X^T Y$を特異値分解します。

「特異値分解」って何だよ！って思われるかもしれません。

これは任意の行列$A$を$U\Sigma V^T$と分解できる事が数学的に証明されています。
特異値分解は上記の分解をする操作です。ここで$U$と$V$はユニタリ行列です。ユニタリ行列についてはググれ。
https://manabitimes.jp/math/2617
また、ここで$\Sigma$は対角成分以外は0の行列です。ただし、正方行列とは限りません。

$SVD(X^T Y) = U \Sigma V$と書くとこうなりますね…
$$Cov(Z_x, Z_y)=w_x^T X^T Y w_y=w_x^T U \Sigma V w_y$$
さて、特異値分解の定義より、$U$と$V$はユニタリ行列です。慣習として$\Sigma$は大きい方から小さい方へと並んでいきます。ここで、仮に$R$という変数を作って

• $R_x = w_x U$
• $R_y = V w_y$

としましょう。すると下記のようになります。
$$w_x^T X^T Y w_y = R_x \Sigma R_y$$
ただし、この時、

• $R_x^T R_x = I$
• $R_y^T R_y = I$

となります。なぜならユニタリ行列と任意の何かを掛け算しても傾きがかわるだけだからです。それは数学の公式になっています。実際、計算すればそうなります。つまり、この最大化問題はこうなります。
$$Cov(Z_x, Z_y) = Rx \Sigma Ry$$
これ…あきらかに$\Sigma$の一番大きい値に1を掛け算するのが最大って思えますよね。つまり$P_x$と$P_y$は一番右上だけが1になる行列です。あとは…$w_x$を求めるだけ。
$$U^T R_x = w_x$$
そうです、$w_x$は$U$のてっぺんだけ取ってきたものなのです!あとは、$w_x$や$w_y$からZは計算できるからいいですね。
そんなわけで、$w_x$と$w_y$は一行だけの行列になりました。

pythonで書くとおおよそこうなります。

    def _fit1(self, x, y):
        u, s, v = np.linalg.svd(x.T @ y, full_matrices=True)
        x_weight, y_weight = (n[:, 0:1] for n in (u, v.T))
        x_score = x @ x_weight
        y_score = y @ y_weight
        self._flip_weight([x_weight, y_weight, x_score, y_score])
        x_load = (x_score.T @ self.x) / (x_score.T @ x_score)
        y_load = (y_score.T @ self.y) / (y_score.T @ y_score)
        self.x -= x_score @ x_load

途中でflipとかしているのは、特異値分解の左右は一意にきまらない事が知られているので一意にする操作をしているのです。

さて、一つだけ問題があります。$Z_y$は1行の行列でした。そんな一部だけ取ってきた行列を解として大丈夫か？実は全然大丈夫じゃないので後述します。

ちゃんとした行列にする

問5. nipalsおよび特異値分解によって求めた部分的最小二乗回帰の潜在変数を全て求めよ。

上記、nipalsと特異値分解の両方のやりかたで$w_x$を1行ずつ作る事に成功しました。ここで特異値分解の場合はさらに2通りに分れます。前提の式を思い出しましょう。

• $X = Z_x X_L + E$
• $Y = Z_y Y_L+ F$

このなかで$E$が無視されていますね！先程の1行だけしか$w_x$が求められず、残りがなかったのは$E$を無視していたからなのです！$E$も$F$も簡単に計算できますよね。
$$E = X - Z_x X_L$$

次からの行列は$E$を$X_2$とおいて計算すればいいのです。ただ、$X_L$の導出ははじめに語った割り算が必要です。

$X = Z_x X_L + E$
これのEが小さいので
$X_L = Z_x X(Z_x^T Z_x)^{-1}$
これで$E$を求める事ができますね!

特異値分解の場合は一寸おまけがあります。
$$A = U\Sigma V^T$$
のうち、$U$と$V^T$の1行目のみ使っていました。もったいなくないですか？折角だから2行目以降もつかおうぜ！という思想があります。解は違うのですが、これは簡易版ですね。実装があります。

過学習

実は部分的最小二乗回帰はかなり性能がいいです。という事は…過学習うるということです。どんな風に過学習するかというと、$E$を延々と限りなく0にしていく事ができるという事です。いけませんね!
なので、それが正しいのかどうかを確認する必要があります。よく行われるのはK fold法です。詳しくはググって下さい。

余裕があれば後日ここに書きます。
また、余裕があれば後日ソースコードも載せます。
ちなみに、ソースコードはscikit-learnとほぼ同じ結果を出す事を確認済みです。

参考にしたもの

おなじみwikipedia。あまり詳しくなかったです。

ここは、ラグランジュの未定乗数法を使った導出をされています。

特異値分解に関して分りやすかったです。

古い論文ですが、nipalsを図解していたのです。良き。

良い子はソースコードを読みましょう。