2014-07-26

次元を上げる

複比配置空間

2次元空間の3個の点の組は1次元射影空間で射影変換的にすべて同値だった
そのときに３点を0,1,無限大に写す行列を考えた
n次元に上げよう

# 適当に作る
k1 <- rnorm(1)
k2 <- rnorm(1)
k3 <- rnorm(1)
x1 <- rnorm(1)
y2 <- rnorm(1)
x3 <- rnorm(1)

P <- matrix(c(k1*x1,k1,k2,k2*y2,k3*x3,k3),2,3)
P
my.st.proj <- function(P){
	M <- P[,1:2]
	M.inv <- solve(M)
	P. <- M.inv %*% P
	P.. <- P.
	N <- diag(c(1/P.[1,3],1))
	L <- diag(c(1,1/P..[2,3]))
	L %*% N %*% M.inv
}
my.st.proj(P)
my.st.proj(P) %*% P
# P はn*(n+1)行列
my.st.proj.n <- function(P){
	n <- length(P[,1])
	M <- diag(rep(1,n-1))
	M <- rbind(M,rep(1,n-1))
	M <- cbind(c(rep(0,n-1),1),M)

	M2 <- P[,1:(n)]
	M2.inv <- solve(M2)
	M.inv <- M %*% M2.inv
	P. <- M.inv %*% P
	P.. <- P.
	N <- diag(c(1/P.[1,n+1],rep(1,n-1)))
	L <- diag(c(rep(1,n-1),1/P..[n,n+1]))
	L %*% N %*% M.inv
}
my.st.proj.n(P)
my.st.proj.n(P) %*% P

n <- 3
Q <- matrix(rnorm(n*(n+1)),n,n+1)
my.st.proj.n(Q)
round(my.st.proj.n(Q)%*%Q,4)

これを使うと、k変数の連立１階常微分方程式の解が射影されてk-1次元で等間隔解析されているときに、k+2個の連続するk-1次元ベクトルのデータには、k-1次元ベクトル保存量があることがわかる
それはk次元空間のk+2個の点をk-1次元射影空間に落とすとk+1個の点は線形変換的に同値で残りの(k+2)-(k+1)個の点が商空間の代表を表す点となっていて、そのk-1次元射影空間の非同次座標が保存されているということ
- k=2のときは、２変数の連立１階常微分方程式の等間隔写像を１次射影直線に作ると、２＋２＝４点について、１次射影空間上の座標(１次元)が保存されていることの拡張版

-Rで確かめておく

my.st.proj.n2 <- function(P){
	n <- min(dim(P))
	P. <- P[1:n,1:n]
	inv.P. <- solve(P.)
	ret <- inv.P.
	#s <- which(P[,n+1]!=0)
	#tmp.s <- s[1]
	tmp.s <- n
	s.others <- (1:n)[-tmp.s]
	tmp.sum <- sum(P[s.others,n+1] * ret[n,s.others])
	ret[n,tmp.s] <- -tmp.sum/P[tmp.s,n+1]
	P.. <- ret %*% P.
	
	for(i in 1:(n-1)){
		ret[i,] <- ret[i,]/P..[i,i]*P..[n,i]
	}
	ret
}
my.st.proj.n2(P)
my.st.proj.n2(P) %*% P

n <- 3
Q <- matrix(rnorm(n*(n+1)),n,n+1)
my.st.proj.n2(Q)
round(my.st.proj.n2(Q)%*%Q,4)

my.matrix.eigen <- function(lambdas,vs){
  Vs <- t(t(vs) * lambdas)
	Vs %*% solve(vs)
}

exp.m <- function(A,n){
  # 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(exp(eigen.out[[1]]*n))
	#B <- diag(eigen.out[[1]]^n)
	X <- V%*%B%*%U
	return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}

my.doubleratio <- function(x){
	if(!is.matrix(x)){
		x <- matrix(x,nrow=length(x))
	}
	P <- x[1:(length(x[,1])-3),]
	Q <- x[4:length(x[,1]),]
	S <- x[2:(length(x[,1])-2),]
	R <- x[3:(length(x[,1])-1),]
	((P-R)/(Q-R))/((P-S)/(Q-S))
}

k <- 3
k2 <-2
A <- matrix(rnorm(k^2),k,k)
A2 <- matrix(rnorm(k2^2),k2,k2)
	# 等間隔なt
t <- seq(from=-1,to=5,length=30)
X <- matrix(0,length(t),length(A[,1]))
X2 <- matrix(0,length(t),length(A2[,1]))
X.init <- runif(k)
X.init2 <- runif(k2)*10
for(i in 1:length(t)){
	X[i,] <- Re(exp.m(A,t[i])[[1]] %*% X.init)
	X2[i,] <- Re(exp.m(A2,t[i])[[1]] %*% X.init2)
}
plot(X)
X.r <- X
plot(X.r)
X. <- t(t(X.r[,1:(k-1)])/X.r[,k])

plot(X.)
#X.. <- matrix(0,length(X.[,1])-k+1,length(X.[1,]))

tmp.list <- list()
for(i in (k+2):(length(X.[,1]))){
	tmp.mat <- t(X.[(i-k-1):(i),])
	tmp.mat.2 <- rbind(tmp.mat,rep(1,length(tmp.mat[1,])))
	tmp <- my.st.proj.n2(tmp.mat.2) %*% tmp.mat.2
	tmp2 <- t(tmp)/tmp[length(tmp[,1]),]
	tmp2[k+1,] <- tmp[,k+1]
	tmp.list[[i-k-1]] <- tmp2
	print(round(tmp2,4))
#my.st.proj.n(t(X.[2:(k+1),])) %*% t(X.[2:(k+1),])
}

Y <- matrix(0,length(tmp.list),k-1)
for(i in 1:length(Y[,1])){
	Y[i,] <- tmp.list[[i]][k+2,1:k-1]
}
plot(Y)

for(i in 1:(length(tmp.list)-1)){
	print(tmp.list[[i]]-tmp.list[[i+1]])
}
par(mfcol=c(1,2))
plot(x,y,type="l")

plot(v,ylim=c(-max(abs(v)),max(abs(v))),type="l")
par(mfcol=c(1,1))


v <- cbind(x,y)
tmp.list <- list()
k <- 4
for(i in k:(length(v[,1]))){
	tmp <- my.st.proj.n(t(v[(1+i-k):(i),])) %*% t(v[(1+i-k):(i),])
	tmp2 <- t(tmp)/tmp[length(tmp[,1]),]
	tmp.list[[i-k+1]] <- tmp2
	print(tmp2)
#my.st.proj.n(t(X.[2:(k+1),])) %*% t(X.[2:(k+1),])
}
for(i in 1:(length(tmp.list)-1)){
	print(tmp.list[[i]]-tmp.list[[i+1]])
}