多次元球面上の点が作る「角度」

  • 任意次元(次元d)の空間に単位球面があり、その上に2<=n<=d個の点があるとする
  • このn個の点が作る「角度」を考える
  • d=n=2の場合は、普通の意味での角度であって、\cos{\theta} = (v1,v2)と言うように内積を取ればよい
  • d=n>2の場合は、n点が作る球面n角形(球面n頂点多面体)の「内側」の球面部分面積の、球面全面積に対する割合を考えればよい
    • それを算出するには、d次元単位球面乱点を発生させ、それをnベクトルの線形和で表したときの、係数が、全部正であるかそうでないかの計算により推定できる
  • では、n
    • n個の単位ベクトルのペアワイズ角度(いわゆる角度)を算出し、それをもとに、d2=n次元空間にn個のベクトルを再配置することができる。ペアワイズ角度は\frac{n(n-1)}{2}個の値を提供し、n次元空間へのベクトル配置は、第1のベクトルの座標を固定すれば、\frac{n(n-1)}{2}個の座標値を逐次確定することであるから、一意に定めることができる
    • d2=nであれば、「角度」は定まる
# Utility functions
# Generate random points on spheres
# df : dimension; r : radius, n : number of points to be generated
RandomSphere<-function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n)
    rs/sqrt(apply(rs^2, 1, sum)) * r
}

# From pairwise cosine matrix X, return k-dimensional coordinates, satisfying X
VectorsFromPairwise<-function(X,k){
	n<-length(X[,1])
	M<-matrix(0,n,n)
	M[1,1]<-1
	M[2,1]<-X[2,1]
	M[2,2]<-sqrt(1-M[2,1]^2)
	for(i in 3:n){
		M[i,1]<-X[i,1]
		for(j in 2:(i-1)){
			M[i,j]<-1/M[j,j]*(X[i,j]-sum(M[i,1:(j-1)]*M[j,1:(j-1)]))
		}
		M[i,i]<-sqrt(1-sum(M[i,]^2))
	}
	# k is enough per argument
	M<-M[,1:k]
	M
}
# n vectors in d (>=n) dimensional space can be placed in n dimensional space 
# with pairwise cosines identical
# vs: n-row x d-column matrix with n vectors
my.vs2sphereSurface <- function(vs){
	d <- length(vs[1,])
	n <- length(vs[,1])
	
	vs <- vs/sqrt(apply(vs^2,1,sum))
	coss <- vs %*% t(vs)
	diag(coss) <- 1
	return(VectorsFromPairwise(coss,n))
}

# Angle is estimated where angle means the area k vectors demarcate on the surface
# vs : k vectors in k dimensional space
# n : number of random points on unit shepre
my.general.angle <- function(vs,n=10^6,ratio=TRUE){
	L <- sqrt(apply(vs^2,1,sum))
	vs <- vs/L
	d <- length(vs[1,])
	R <- RandomSphere(df=d,n=n)
	# To calculate linear expression coefficients for random points with vs
	M <- solve(t(vs))
	K <- M %*% t(R)
	tmp <- sum(apply(K>0,2,prod))/n
	if(ratio){
		return(tmp)
	}else{
		S <- 2*pi^(d/2) / gamma(d/2) # Area of sphere in d-dimensional space
		return(tmp*S)
	}
}

# example
# Generate random vectors
k<-5 # dimension
n<-30 # number of unit vectors
X<-RandomSphere(df=k,n=n)
X
# pairwise cosine matrix

Y<-X%*%t(X)
# diagonal values should be 1
diag(Y)<-1
Y

# calculate coordinate
M<-VectorsFromPairwise(Y,k)
M
# Pairwise cosines are identical before and after the coordinate transformation
Y2<-M%*%t(M)
max(abs(Y-Y2))


# d dimension, n vectors
d <- 10
n <- 5
vs <- RandomSphere(df=d,n=n)

# pairwise cosines
vs %*% t(vs)

# new coordinates in n dimensional space
new.vs <- my.vs2sphereSurface(vs)

# pairwise cosines should be the same as before
new.vs %*% t(new.vs)

my.general.angle(new.vs)

# In case of 2 dimension, angle estimation works fine
theta <- pi/2
theta2 <- pi/6

vs <- rbind(c(cos(theta),sin(theta)),c(cos(theta+theta2),sin(theta+theta2)))

my.general.angle(vs)

# dim = 3...
vs <- rbind(c(1,0,0),c(0,1,0),c(0,0,1))

my.general.angle(vs)