2013-12-03

周期性と巡回行列

周期性というのは時間に関して、時刻間の差によって観察値の違い(ばらつき)が規定され、その時刻の差は周期のmoduloとして定められる者と言える
時刻ごとの観察を変数として並べ、それらの関係を分散共分散行列として表すと、その行列は巡回行列になる
- Rで作ればこんな感じ
  - まずは、どのセルにどんなのが入っているかに注意して作る

# 変数の数
n <- 6
L <- round(runif(n),2)
sigma <- matrix(0,n,n)
for(i in 1:n){
	tmp <- (1:n+i-2)%%n+1
	sigma[cbind(tmp,1:n)] <- L[i]
}
L
sigma

> L
[1] 0.10 0.56 0.30 0.68 0.71 0.01
> sigma
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.10 0.01 0.71 0.68 0.30 0.56
[2,] 0.56 0.10 0.01 0.71 0.68 0.30
[3,] 0.30 0.56 0.10 0.01 0.71 0.68
[4,] 0.68 0.30 0.56 0.10 0.01 0.71
[5,] 0.71 0.68 0.30 0.56 0.10 0.01
[6,] 0.01 0.71 0.68 0.30 0.56 0.10

- 次に巡回行列をこんなルールで作る
  - - ただし $P$ は $(i,i-1)$ と $(1,n)$ に１が立ち、それ以外は0の巡回行列

P <- diag(rep(1,n))
P <- P[,c(2:n,1)]
P
sigma2 <- diag(L[1],n)
tmpP <- diag(rep(1,n))
for(i in 2:n){
	tmpP <- tmpP %*% P
	sigma2 <- sigma2 + L[i]*tmpP
}
sigma2
sigma

> P
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    0    0    1
[2,]    1    0    0    0    0    0
[3,]    0    1    0    0    0    0
[4,]    0    0    1    0    0    0
[5,]    0    0    0    1    0    0
[6,]    0    0    0    0    1    0
> sigma2
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.10 0.01 0.71 0.68 0.30 0.56
[2,] 0.56 0.10 0.01 0.71 0.68 0.30
[3,] 0.30 0.56 0.10 0.01 0.71 0.68
[4,] 0.68 0.30 0.56 0.10 0.01 0.71
[5,] 0.71 0.68 0.30 0.56 0.10 0.01
[6,] 0.01 0.71 0.68 0.30 0.56 0.10
> sigma
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.10 0.01 0.71 0.68 0.30 0.56
[2,] 0.56 0.10 0.01 0.71 0.68 0.30
[3,] 0.30 0.56 0.10 0.01 0.71 0.68
[4,] 0.68 0.30 0.56 0.10 0.01 0.71
[5,] 0.71 0.68 0.30 0.56 0.10 0.01
[6,] 0.01 0.71 0.68 0.30 0.56 0.10

ここで上でLとした、巡回行列のセルの値を $c=(c_0,c_1,...,c_{n-1})=(L[1$ ,L[n],L[n-1],...,L[2])]として
これの固有値分解をしたものは、のn個の根を使ってで得られるという
- Rでやってみる

eigen.out <- eigen(sigma)
ws <- exp(2*pi*1i*(0:(n-1))/n)
lambdas <- rep(0,n)
for(i in 1:n){
	lambdas[i] <- sum(L[c(1,(n):2)] * ws[i]^(0:(n-1)))
}
eigen.out[[1]]
ord <- order(Mod(lambdas),decreasing=TRUE)
lambdas[ord]

> eigen.out[[1]]
[1]  2.36+0.0000000i -0.01+0.8313844i -0.01-0.8313844i -0.80+0.1212436i
[5] -0.80-0.1212436i -0.14+0.0000000i
> ord <- order(Mod(lambdas),decreasing=TRUE)
> lambdas[ord]
[1]  2.36+0.0000000i -0.01+0.8313844i -0.01-0.8313844i -0.80-0.1212436i
[5] -0.80+0.1212436i -0.14+0.0000000i

こうしてみてくると巡回行列は、ぐるり１順を表す行列Pのべきに係数重みをつけた行列であって、"unitary discrete Fourier transform matrix"というのを使うと、任意の巡回行列が、そのunitary discrete Fourier transform matrix (Un)とそのユニタリ行列とそれに挟まれた対角行列の積に分解できることがわかり(これは固有値分解)、その対角行列もfourier transform matrixと元の順が居行列の第１列(それに巡回行列の情報はすべて載っている)との積で簡単に得られる、という構造になっている
- Rでやってみる

# 正方行列を作る。fourier transform matrixやその部品
U <- F <- matrix(0,n,n)
U.address <- which(U==0,arr.ind=TRUE)-1
# この部分の計算が、「２つの時刻j,jを用いたフーリエ的計算」になっている
f.jk <- exp(-2*U.address[,1]*U.address[,2]*pi*1i/n)
F <- matrix(f.jk,n,n)
U <- F/sqrt(n)
# 巡回行列sigmaの第１列だけ使っている
S <- F %*% matrix(sigma[,1],ncol=1)
# 固有値分解を再計算して、元の巡回行列と比較、計算機誤差以外は一致
t(Conj(U)) %*% diag(c(S)) %*% U
sigma

> sigma
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.10 0.01 0.71 0.68 0.30 0.56
[2,] 0.56 0.10 0.01 0.71 0.68 0.30
[3,] 0.30 0.56 0.10 0.01 0.71 0.68
[4,] 0.68 0.30 0.56 0.10 0.01 0.71
[5,] 0.71 0.68 0.30 0.56 0.10 0.01
[6,] 0.01 0.71 0.68 0.30 0.56 0.10
> t(Conj(U)) %*% diag(c(S)) %*% U
        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
[1,] 0.10-0i 0.01-0i 0.71+0i 0.68+0i 0.30-0i 0.56+0i
[2,] 0.56+0i 0.10-0i 0.01-0i 0.71+0i 0.68+0i 0.30-0i
[3,] 0.30+0i 0.56+0i 0.10-0i 0.01-0i 0.71+0i 0.68+0i
[4,] 0.68-0i 0.30-0i 0.56+0i 0.10-0i 0.01-0i 0.71+0i
[5,] 0.71+0i 0.68+0i 0.30+0i 0.56-0i 0.10-0i 0.01+0i
[6,] 0.01-0i 0.71+0i 0.68+0i 0.30+0i 0.56+0i 0.10-0i