ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2016-12-27 統計遺伝学分野で生き抜くための数学〜精選7〜

[][][][][]行列で最善の解

---
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)
```


# Exercise 1

ムーアペンローズ疑似逆行列の解の幾何的な意味を以下の手順で確認せよ。

## Exercise 1-1

$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$のとき
の2直線を描き、その交点をMASS::ginv()関数を求め、その点$(\hat{a_1},\hat{a_2})$が交点にあることを、打点することによって示せ。

## Exercise 1-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}
$$

これは3直線の場合である。

$\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}$とする。


3直線を描図し、ムーアペンローズ疑似逆行列による解$(\hat{a_1},\hat{a_2})$を打点せよ。

## Exercise 1-3

Excerise 1-3 の例で、第1の直線は$(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本の直線を描け。

同様に、第1、第2、第3の直線と、それに対応する、$(\hat{a_1},\hat{a_2})$を通る3直線を描け。

## Exercise 1-4

第1の直線は、$(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})$である。
この垂直な直線が第1の直線と交わる交点$P_1$を求めよ。

$(\hat{a_1},\hat{a_2})$と交点$P_1$との距離を求めよ。

同様に、行列が定める3つの直線のそれぞれについて、行え。

$||y_i,\hat{y_i}||$とこの距離の関係は何か。

## Exercise 1-5

$(a_1,a_2)$平面の任意の点に関して、その点から、直線への垂線の足を求める関数を作成せよ。

その関数を用いて、行列が定める3直線への足、3点が求まる。
$||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})$の意味を確認せよ。

## Exercise 1-6

$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本引ける。
ムーアペンローズ疑似逆行列の解は直線上の1点である。

$||(a_1,a_2)||$が原点から最短距離になっていることを確かめよ。

## Exercise 1-7

$n=m$,$n>m$,$n<m$の3つの場合について、ムーアペンローズ疑似逆行列の解の幾何学的な意味を説明せよ。

## Exercise 1-8

$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$でコントロールする。

ムーアペンローズによる解の選び方と、正則化手法での解の選び方との類似性についてコメントせよ。