---
title: "Best answer with Matrix: Moore-Penrose Pseudo-inverse 行列で最善の解: ムーアペンローズ疑似逆行列"
author: "ryamada"
date: "2016年12月27日"
output:
html_document:
toc: true
toc_depth: 6
number_section: true
---
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```
変数の数と等式の数が一致しているとき、うまく逆行列が取れれば
$$
y = X \mathbf{a}
$$
はきっちりと解けて$\mathbf{a}$が一意に求まる。
その際、
$$
\mathbf{a} = X^{-1} y = (X^TX)^{-1}X^T y
$$
であった。
線形回帰の場合は、サンプルの数(ベクトル$y$の長さ:$n$とする)に対して、説明変数の数$X$の列数($m$とする)が小さいので
$$
y \sim X \mathbf{a}
$$
として、
$$
\sum_{i=1}^n (y_i-\hat{y_i})^2
$$
が最小になるような$\mathbf{a}$を推定するわけだが、そのとき
$$
\mathbf{a} = (X^T X)^{-1} X^{-1}
$$
で求まるのだった。
結局、行列$X$の列数によらず、$\mathbf{a}$について、「これが一番」という答えが
$$
\mathbf{a} = (X^TX)^{-1}X^T y
$$
で求まる。
ここで、「これが一番」というのは、最小二乗が0であったり、「かっちり連立方程式が解ける」ということだったりした。
行列$X$の行数$n$,列数$m$について
+ $n=m$の場合：連立方程式を解く
+ $n>m$の場合：線形回帰をする
だった。
網羅されていないのは
+ $n<m$の場合である
この場合にも、うまく行く方法があり、これをムーアペンローズ疑似逆行列と読んだり、一般化逆行列と読んだり、疑似逆行列と呼んだりする。
$n < m$の場合には、$(X^TX)^{-1}$がうまく計算されないので、別の式を用いる。
Xの行数が列数より小さい場合には、転置して考える。
転置($X'=X^T$)すると、そのムーアペンローズ疑似逆行列は、$(X'^TX')^{-1}X'^T$と計算できるので、それを再度転置すれとうまく行く($((X'^TX')^{-1}X'^T)^T$)。
結局、覚えるべきは
+ $n \ge m$ のときは $gen.inv(X)=(X^TX)^{-1}X^T$
+ $n < m$のときは、$X_{gen.inv} = (gen.inv(X^T))^T$
$n < m$の場合を$X^T(X^TX)^{-1}$と書くこともある。
```{r}
n <- 3
m <- 2 # m < n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
```
```{r}
n <- 3
m <- 3 # m = n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
```
```{r}
n <- 2
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
# solve(t(X)%*%X) %*% t(X) # Error
tX <- t(X)
t(solve(t(tX)%*%tX)%*%t(tX))
t(X)%*%solve(X%*%t(X))
```
Rでは、MASS パッケージにginv() 関数があり、Xの行数・列数に関係なく、ムーアペンローズ疑似逆行列を算出できる。
任意の行列を次のように分解することを特異値分解と言うが
$$
X = U\Sigma V^T
$$
ムーアペンローズ疑似逆行列は
$$
gen.inv(X) = V \Sigma^{+} U^T
$$
で与えられる。ただし$\Sigma$は対角成分に特異値を持つ行列であり、$/Simga^{+}$は特異値の逆数を対角成分に持つ。
実際、Rではこれを利用して、MASSパッケージのginv()関数(g:generalized, inv:inverse)がXの行数・列数を場合分けすることなく、ムーアペンローズ疑似逆行列を返す。
```{r}
library(MASS)
n <- 3
m <- 2 # m < n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
ginv(X)
```
```{r}
n <- 3
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
ginv(X)
```
```{r}
n <- 2
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
tX <- t(X)
t(solve(t(tX)%*%tX)%*%t(tX))
ginv(X)
```
ムーアペンローズ疑似逆行列の解の幾何的な意味を以下の手順で確認せよ。
$n=m$
$$
\begin{pmatrix}y_1\\y_2\end{pmatrix} = \begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\end{pmatrix}\begin{pmatrix}a_1\\a_2\end{pmatrix}
$$
今、$a_1,a_2$の値を求めるというのは、$(a_1,a_2)$を次元座標として、2次元平面にある2本の直線の交点座標を求めることである。
$$
x_{11} a_1 + x_{12}a_2=y_1\\
x_{21} a_1 + x_{22}a_2=y_2
$$
今、$y_1=3,y_2=1,x_{11}=4,x_{12}=2,x_{21}=-1,x_{22}=3$のとき
の２直線を描き、その交点をMASS::ginv()関数を求め、その点$(\hat{a_1},\hat{a_2})$が交点にあることを、打点することによって示せ。
$n > m$
$$
\begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} = \begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\\x_{31},x_{32}\end{pmatrix}\begin{pmatrix}a_1\\a_2\\end{pmatrix}
$$
これは３直線の場合である。
$\begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} = \begin{pmatrix}1\\2\\3\end{pmatrix}$,$\begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\\x_{31},x_{32}\end{pmatrix} = \begin{pmatrix}2,3\\1,-2\\-3,-4\end{pmatrix}$とする。
３直線を描図し、ムーアペンローズ疑似逆行列による解$(\hat{a_1},\hat{a_2})$を打点せよ。
Excerise 1-3 の例で、第１の直線は$(a_1,a_2)$座標について、$y_1 = x_{11}a_1 + x_{12}a_2$で表されているのに対し、解$(\hat{a_1},\hat{a_2})$は、別の直線、$\hat{y_1} = x_{11}\hat{a_1} + x_{12}\hat{a_2}$を表している。
2本の直線を描け。
同様に、第１、第２、第３の直線と、それに対応する、$(\hat{a_1},\hat{a_2})$を通る３直線を描け。
第１の直線は、$(a_1 = \frac{y_1}{x_{11}},0)$と$(0,\frac{y_1}{x_{12}})$を通る直線である。
それに対応する$(\hat{a_1},\hat{a_2})$を通る直線は$(a_1 = \frac{\hat{y_1}}{x_{11}}=\frac{\hat{a_1}}{x_{11}},0)$と$(0,\frac{\hat{y_1}}{x_{12}}=\frac{\hat{a_2}}{x_{12}})$を通る直線である。
これらは平行である。
点$(\hat{a_1},\hat{a_2})$を通り、この2直線と垂直な直線を引け。その直線の傾きは、$(x_{12},-x_{11})$である。
この垂直な直線が第１の直線と交わる交点$P_1$を求めよ。
$(\hat{a_1},\hat{a_2})$と交点$P_1$との距離を求めよ。
同様に、行列が定める３つの直線のそれぞれについて、行え。
$||y_i,\hat{y_i}||$とこの距離の関係は何か。
$(a_1,a_2)$平面の任意の点に関して、その点から、直線への垂線の足を求める関数を作成せよ。
その関数を用いて、行列が定める３直線への足、３点が求まる。
$||y_i,\hat{y_i}||$がわかるので、$\sum_{i=1}^3 (y_i-\hat{y_i})^2$も求まるはずである。
この値が$(a_1,a_2)$平面に、どのような高低を作るかを図示し、$(\hat{a_1},\hat{a_2})$の意味を確認せよ。
$n < m$
$$
\begin{pmatrix}y_1\end{pmatrix} = \begin{pmatrix}x_{11},x_{12} \end{pmatrix}\begin{pmatrix}a_1\\a_2\end{pmatrix}
$$
この場合、直線は１本引ける。
ムーアペンローズ疑似逆行列の解は直線上の1点である。
$||(a_1,a_2)||$が原点から最短距離になっていることを確かめよ。
$n=m$,$n>m$,$n<m$の３つの場合について、ムーアペンローズ疑似逆行列の解の幾何学的な意味を説明せよ。
$n<m$の場合に、$y=X\mathbf{a}$の解は、直線上のいずれの点でもよかったが、ムーアペンローズ疑似逆行列は、原点からの距離が最短になるような点$||(\hat{a_1},\hat{a_2})||$を選んだ。
正規化手法と呼ばれる手法では、
$$
\sum_{i=1}^n (y_i-\hat{y_i})^2 + \lambda ||\mathbf{a}||^p
$$
を最小にすることで、残差$\sum_{i=1}^n (y_i-\hat{y_i})^2$を小さくことのみを目指すのではなく、$||\mathbf{a}||^p$という形で、解の取り方に制約を持たせること、その制約の強さをパラメタ$\lambda$でコントロールする。
ムーアペンローズによる解の選び方と、正則化手法での解の選び方との類似性についてコメントせよ。
