---
title: "Properties of matrices 行列の特徴"
author: "ryamada"
date: "2016年12月22日"
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)
```
When we use matrices for statistical analyses, we change bases so that "better" axes should grab data sets in simpler or better ways.
The change of bases are achieved with linear transformation.
Therefore matrices that are transformed each other by linear transformation should be considered "similar".
Among the similar matrices, the simplest one tends to represent data structure the "best" way from data interpretation standpoint.
The simplest one among the similar matrices is Jordan normal.
Or the simplest one is diagonal, if diagonalizable, and, in many cases the matrices in statistical analyses are diagonizable.
統計解析に行列を使うとき、データセットが単純に理解できるような基底の取り直しをすることがよくある。
基底の取り直しは線形変換で行うので、線形変換で互いに行きあえる行列同市は「相似(似ている)」関係にあると言う。
相互に相似な関係にある行列の中で、最も単純な行列を代表とすることがあり、それを使うことでデータ解釈が最もわかりやすくなる。
相似な行列の中で最も単純なのは、ジョルダン標準形のそれである。
対角化可能な行列なら、ジョルダン標準形は対角化行列である。実際、統計解析で現れる行列の多くは、対角化可能である。
任意の正方行列$M$はジョルダン標準形$J$と逆行列を持つ行列$R$とを用いて、次のように表せる。
$$
M = R J R^{-1}
$$
$M' = R'JR'^{-1}$となるとき、$M$と$M'$は相似であると言う。
以下では、ジョルダン標準形が対角行列であるようなもの(対角化可能行列)に絞って話を進めるが、特に断らなければ、提示している性質等は対角化可能行列に限るものではない。
相互に相似な行列は見かけは違うが、「本質的に同じ」なので、その「本質」に関する性質は共通である。
+ 階数(Rank)
+ 行列式(Determinant)
+ トレース(跡) (Trace)
+ 固有値(Eigen values)
Rank, Determinant, Trace, Eigen valuesは次のように計算できる
```{r}
library(Matrix)
d <- 4
M <- matrix(rnorm(d^2),d,d)
eigen(M)[[1]]
rankMatrix(M)[[1]]
det(M)
sum(diag(M)) # Trace
```
逆行列を持つ行列Rを使って、対角行列Jを変換した行列M($M=RJR^{-1}$)がJと相似であることを、固有値が一致していることを確認して確かめよ。
そして、ランク、行列式、トレースが同一であることを確かめよ。
```{r}
d <- 4
J <- diag(rnorm(d))
R <- matrix(rnorm(d^2),d,d)
M <- R %*% J %*% solve(R)
eigen(J)[[1]]
eigen(M)[[1]]
rankMatrix(J)[[1]]
rankMatrix(M)[[1]]
det(J)
det(M)
sum(diag(J)) # Trace
sum(diag(M)) # Trace
```
行列が相互に相似かどうかは、(すべての固有値が異なるならば)固有値のセットが等しいことで確認できる。
したがって相似な行列に共通する特徴は、固有値と紐づけて算出できるはずである。
以下を確かめよ。
$$
determinant = \prod_{i=1}^n \lambda_i\\
trace = \sum_{i=1}^n \lambda_i
$$
$d$個の異なる固有値を持つ行列を作り、それと相似な行列を複数作れ。
それらの固有値のセットが一致すること、トレース・行列式が一致することを確かめよ。
また、トレースの計算は、対角成分の和を取る方法と、固有値の和を取る方法の両方を実施せよ。
行列式の計算もdet()関数を使用する方法と固有値の積を取る方法の両方を実施せよ。
```{r echo=FALSE}
lambda <- rnorm(d)
M <- diag(lambda)
M.list <- list()
n <- 10
for(i in 1:n){
R <- matrix(rnorm(d^2),d,d)
M.list[[i]] <- R %*% M %*% solve(R)
}
lapply(M.list,det)
lapply(M.list,diag)
lapply(M.list,function(M){eigen(M)[[1]]})
```
$M=M^T$である正方行列を対称行列と言う。
対称行列は
$$
M = R J R^{-1} = R J R^T\\
R^T=R^{-1}\\
R R^T = R^T R = I
$$
と対角化で現れる行列$R$が特別な性質(直交行列)を持つ。
確かめてみる。
```{r}
d <- 5
n.iter <- 10
for(i in 1:n.iter){
M <- matrix(rnorm(d^2),d,d)
M <- M + t(M) # 対称行列化
eigen.out <- eigen(M)
print(round(eigen.out[[2]] %*% t(eigen.out[[2]]),8))
}
```
対角行列と直交行列によって相似な関係にある行列を複数作り、それらが対称行列であることを確かめよ。
なお、直交行列をランダムに発生するには、GPRrotation::Rnadom.Start()関数を用いるとよい。
```{r echo=FALSE}
library(GPArotation)
d<-5
n.iter <- 10
for(i in 1:n.iter){
lambdas <- rnorm(d)
J <- diag(lambdas)
#s <- sample(1:(d-1),sample(1:(d-1),1))
#for(i in s){
# J[i,i+1] <- 1
#}
R <- Random.Start(d)
print("unit?")
print(round(R %*% t(R),8))
M <- R %*% J %*% t(R)
print("zero?")
print(round(M - t(M),8)) # 対称行列なら0行列
}
```
行列には相似という間柄があり、それは固有値で決まることが分かった。
また、相似な行列の間で保存される性質があることもわかった。
また、$R M R^{-1}$という変換が相似な行列を作ることもわかった。
また対象行列の場合 $R M R^{-1} = R M R^T$と書け、Rが直交行列であることもわかった。
これらは、一般的な行列の話であったが、この節では、統計解析でよく登場する特別な行列を取り上げ、その場合に固有値・相似行列の特徴がどのような性質を持つかを確認する。
$n$個の要素があって、それらの要素ペアの距離を行列とすることがある。
要素間の関連を調べるときの基本的なデータ構造であり
距離の場合、自身から自身への距離は0である。
自身以外への距離は正である。
AからBへの距離とBからAへの距離は同じなので、距離行列は、次の性質を持つ。
+ 対角成分が0
+ 対称
+ 非対角成分は正
このような行列では、以下のような特徴がある。
```{r}
d <- 5
M <- matrix(abs(rnorm(d^2)),d,d)
diag(M) <- 0
M <- M+t(M) # 対称にする
eigen.out <- eigen(M)
trace <- sum(diag(M))
trace2 <- sum(eigen.out[[1]])
determinant <- det(M)
determinant2 <- prod(eigen.out[[1]])
eigen.out[[1]]
round(eigen.out[[2]] %*% t(eigen.out[[2]]),6)
```
$k$個の変数、$n$個のサンプルがあるとき、$k\times k$の分散共分散行列が算出できる。
対角成分は$k$個の変数の分散である。
トレースは$k$個の変数の分散の和である。
$k$個の変数の分散の和が固有値の和である。
対称行列である。これは距離行列の場合と同じ。
```{r}
k <- 5
n <- 100
X <- matrix(0,n,k)
X[,1] <- rnorm(n)
for(i in 2:k){
X[,i] <- rnorm(n) + X[,i-1]
}
M <- cov(X)
M
round(M-t(M),8)
det(M)
sum(diag(M))
eigen.out <- eigen(M)
round(eigen.out[[2]] %*% t(eigen.out[[2]]),8)
eigen.out[[1]]
```
分散共分散行列は、対象行列であって、対角成分はすべて正。
また、以下に示すように、固有値もすべて正。
```{r}
n.iter <- 1000
im.part <- matrix(0,n.iter,2)
re.part <- matrix(0,n.iter,2)
for(i in 1:n.iter){
d <- sample(2:50,1)
X <- matrix(0,n,k)
X[,1] <- rnorm(n)
for(j in 2:k){
X[,j] <- rnorm(n) + X[,j-1]
}
M <- cov(X)
M
eigen.out <- eigen(M)
im.part[i,] <- range(Im(eigen.out[[1]]))
re.part[i,] <- range(Re(eigen.out[[1]]))
}
matplot(im.part,type="l")
matplot(re.part,type="l",ylim=c(-1,max(re.part)+1))
plot(re.part[,1],ylim=c(-1,max(re.part[,1]+1)))
abline(h=0,col=2)
min(re.part[,1])
```
固有値がすべて正であるような対称行列は正定値行列を呼ばれている。
正定値行列$M$は、任意の0でないベクトル$v$に対して
$$
v^T M v >0
$$
となる。
正定値行列の特殊例でもある、単位行列の場合は、以下のように、「ベクトルの長さの二乗」である。
これが0ベクトルなら0となり、そうでなければ正になる、というのは、ベクトルの「長さ」を定義するのに適切な性質を持つ行列を正定値行列と言う、というようにも理解できる。
逆に言えば、統計解析で登場する行列は、この性質を満たすものが多いことともなる。
$$
v^T M v = v^T v
$$
このようなMは、ベクトルの長さの定め方でもあり、$v^T M v$の値を決めるという意味では、「計量」を定めるとも言える。
対称行列の固有値分解 $M = RJR^{-1}=RJR^T$で、Vは回転行列の性質を持つから、
対角成分がすべて正である対角行列と回転行列から作る行列は正定値行列であることになる。
```{r}
library(GPArotation)
d <- 5
#V <- Random.Start(d)
V <- matrix(rnorm(d^2),d,d)
lambdas <- abs(rnorm(d))
M <- V %*% diag(lambdas) %*% solve(V)
eigen(M)[[1]]
lambdas
```
```{r}
n.vec <- 1000
vecs <- matrix(rnorm(n.vec*d),nrow=d)
range(diag(t(vecs) %*% M %*% vecs))
```
多変量正規分布は、統計解析でよく使われる。
$$
N(\mu,\Sigma)
$$
とパラメタ表現される。
変数の数$n$個の多変量正規分布の場合、
$\mu$は長さ$n$のベクトルであり、各変数の平均値を表す。
$\Sigma$は$n \times n$行列であり、正定値行列であることになっている。
正定値行列は、対角化可能で、$\Sigma = RJR^{-1}$と対角化したときに、$R$が直交行列($R^{-1}=R^T$)であるようなものであった。
これを利用して、多変量正規分布からの乱数発生をしてみる。
```{r}
library(mvtnorm)
#help(rmvnorm)
sigma <- matrix(c(4,2,2,3), ncol=2)
x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
colMeans(x)
var(x)
plot(x)
```
2変数標準正規分布は、３変数の平均がともに0で、分散共分散行列が$2\times 2$の単位行列であるようなもののことである。
```{r}
d <- 2
m <- rep(0,d)
sigma <- diag(rep(1,d))
m
sigma
n <- 10000
x <- rmvnorm(n=n,mean=m,sigma=sigma)
apply(x,2,mean)
apply(x,2,var)
cov(x)
round(cov(x),8)
plot(x,asp=TRUE,pch=20,cex=0.1)
```
2変量正規分布において、各変数の分散を(1,1)からずらして、乱点を発生せよ。
なお、分散共分散行列は対角行列のままとする。
```{r echo=FALSE}
sigma2 <- diag(runif(d))
x <- rmvnorm(n=n,mean=m,sigma=sigma2)
apply(x,2,mean)
apply(x,2,var)
round(cov(x),8)
sigma2
plot(x,asp=TRUE,pch=20,cex=0.1)
```
4-1で用いたSigmaと相似な行列で対称なものは正定値行列であるから、mvtnorm::rmvnorm()のsigma引数として適当なはずである。
そのようなsigmaを作り、乱点を発生せよ。
```{r}
R <- Random.Start(d)
sigma3 <- R %*% sigma2 %*% t(R)
x <- rmvnorm(n=n,mean=m,sigma=sigma3)
apply(x,2,mean)
apply(x,2,var)
round(cov(x),8)
sigma3
plot(x,asp=TRUE,pch=20,cex=0.1)
```
分散共分散行列を $R J R^{-1}$にて変換すると、乱点分布が回転し、楕円の軸が斜めになった。
実際、分散共分散行列の変換に用いた直交行列$R$を利用して、乱点を変換すると、楕円の軸を正すことができる。
以下のコードを読み、変換した分散共分散行列に従う２変量正規分布の軸が補正されていることを、コードの各行の処理にコメントを加えることで確認せよ。
```{r}
R <- Random.Start(d)
sigma3 <- R %*% sigma2 %*% t(R)
x <- rmvnorm(n=n,mean=m,sigma=sigma3)
apply(x,2,mean)
apply(x,2,var)
cov(x)
round(cov(x),8)
plot(x,asp=TRUE,pch=20,cex=0.1)
x2 <- t(t(R) %*% t(x))
cov(x2)
plot(x2,asp=TRUE,pch=20,cex=0.1)
```
次のようにして生成される、多変量の観察データ x があったとする。
```{r}
d <- 3
n <- 1000
x <- matrix(rnorm(d*n),ncol=d)
x <- apply(x,2,cumsum)
plot(as.data.frame(x),pch=20,cex=0.1)
```
```{r rgl=TRUE}
plot3d(x)
```
xの各変数の標本平均を求めよ。
xの分散共分散行列を求めよ。
```{r echo=FALSE}
apply(x,2,mean)
cov(x)
```
4-1 で求めた平均値ベクトルと分散共分散行列とに基づく、多変量正規乱数を発生し、ペアワイズにプロットせよ。
```{r echo=FALSE, rgl=TRUE}
n <- 1000
x.new <- rmvnorm(n=n,mean=apply(x,2,mean),sigma=cov(x))
plot3d(x.new)
```
```{r echo=FALSE, rgl=TRUE}
plot3d(x.new)
points3d(x,col=2,radius=0.1)
```
回転行列を用いて、xを回転し(軸を取り直し)、その標本平均と分散共分散行列とを求め、4-1のそれと比較せよ。
各変数の平均も分散も異なることを確かめよ。
```{r echo=FALSE}
R <- Random.Start(d)
x2 <- t(R %*% t(x))
print("mean x")
apply(x,2,mean)
print("mean x rotated")
apply(x2,2,mean)
print("cov x")
cov(x)
print("cov x rotated")
cov(x2)
```
分散は異なるが、分散の和は変わっていないことを、分散共分散行列のトレースを計算することで確かめよ。
同様に、分散共分散行列の固有値の和としてトレースを計算することでも確かめよ。
```{r echo=FALSE}
sum(diag(cov(x)))
sum(diag(cov(x2)))
sum(eigen(cov(x))[[1]])
sum(eigen(cov(x2))[[1]])
```
データセットxの軸の取り方を回転することによって変えても、各軸の分散の和は変わらない。
分散の和は変わらないが、分散の総和の分配具合は変わる。
そのことを、多数の回転行列をランダムに発生させ、データセットを回転した上で各軸方向の分散を計算し、その分散の値を大小順序でソートしてプロットし、分配具合が変わることを具体的に確認せよ。
```{r echo=FALSE}
n.iter <- 100
sorted.vars <- matrix(0,n.iter,d)
for(i in 1:n.iter){
R <- Random.Start(d)
x.tmp <- t(R %*% t(x))
cov.tmp <- cov(x.tmp)
sorted.vars[i,] <- sort(diag(cov.tmp))
}
matplot(t(sorted.vars),type="l",xlab="order",ylab="var")
```
xの分散共分散行列の固有値を求めることで、分散共分散行列を対角化するような回転をしたときの、各軸の分散を求めることができる。
その分散を4-4のプロットに重ねてプロットせよ。
```{r,echo=FALSE}
x.cov <- cov(x)
lambdas <- sort(eigen(x.cov)[[1]])
matplot(t(sorted.vars),type="l",xlab="order",ylab="var",col=1)
points(1:d,lambdas,col=2,pch=20)
for(i in 1:(d-1)){
segments(i,lambdas[i],i+1,lambdas[i+1],col=2)
}
```
4-4,4-5と同様のプロットを変数の数を増やして実施せよ。
```{r}
d <- 10
n <- 1000
x <- matrix(rnorm(d*n),ncol=d)
x <- apply(x,2,cumsum)
plot(as.data.frame(x[,1:5]),pch=20,cex=0.1)
```
```{r echo=FALSE}
n.iter <- 100
sorted.vars <- matrix(0,n.iter,d)
for(i in 1:n.iter){
R <- Random.Start(d)
x.tmp <- t(R %*% t(x))
cov.tmp <- cov(x.tmp)
sorted.vars[i,] <- sort(diag(cov.tmp))
}
matplot(t(sorted.vars),type="l",xlab="order",ylab="var")
```
```{r,echo=FALSE}
x.cov <- cov(x)
lambdas <- sort(eigen(x.cov)[[1]])
matplot(t(sorted.vars),type="l",xlab="order",ylab="var",col=1,ylim=range(c(sorted.vars,lambdas)))
points(1:d,lambdas,col=2,pch=20)
for(i in 1:(d-1)){
segments(i,lambdas[i],i+1,lambdas[i+1],col=2)
}
```
