多次元で垂線を下ろす - ryamadaの遺伝学・遺伝統計学メモ

N次元空間を考える
ある点Aがある
いくつかの点の集合Xsがある
原点とXsの各点を結んだ線によって「籠」のようになった範囲は錐
Aの方向がこの錐の内部かどうかの判定を昨日した
Xsの頂点を結ぶ平面を考える
頂点に囲まれた内部とその外部を区別する
この平面にAから垂線を下ろす。足をA'とする
Aと平面の距離はこの垂線の長さ
垂線の足A'がXsの頂点内部にあるとき、AとXsが表すN次元空間中のN-1次元線形多様体との距離はこれである
垂線の足A'がXsの外部のときは、まだ、AとXsの囲む多様体との距離はわからない
多様体の「周辺」を作るN-2次元線形多様体が複数あるので、それへの距離を測定する
N次元N点からスタートすれば、N通りのN-1点の組み合わせ、次いで、N-1通りのN-2点の組み合わせと、次元を下げていき、いつか、垂線の足が「内部」に落ちることによって「最短距離」が定まる
その基本ための、「内部判定」が昨日
今日は、「垂線の足」「垂線の長さ」の算出
$\sum_{i=1}^k a_i v_i =p$ というベクトルpがある。pはk個の線形独立なベクトルの線形和
k個のベクトルの重心をCとして、 $v'_i=v_i-C$ とすれば、 $v'_i$ とpとが垂直なとき、それが、pからのこの線形和で表された単体からの垂線
原点から、原点を含まない単体 $\sum_{i=1}^k a_i v_i =p,\sum a_i =1, a_i \ge 0$ への垂線の求め方は、この条件の連立方程式を解くことである。

nrm<-function(A){
	sqrt(sum(A^2))
}

#k次元空間中のk個のベクトルの作る凸包への垂線とその長さを求める

MinDistance2<-function(A,Xs,bothside=TRUE){
 df<-length(A)
 N<-length(Xs[,1])
 Xs1<-t(t(Xs)-A)
 C<-apply(Xs1,2,sum)/N
 Xs2<-t(t(Xs1)-C)
# \sum_{i=1}^N a_i v_i と表したときに \sum_{i=1}^N a_i=1という制約がある
# v_iは与えられている凸包の頂点座標(ただし、重心を原点にとりなおしている
# 垂線の足をbとすると
# Ma=bという連立方程式を解いてaを求めるたい
# ただし凸包はN次元ではなくてN-1次元なので、N次元のうち１次元分は取り除
いて
# その代わり係数の和が１という制約等式を入れれば
# あとは行列を使って連立方程式を解くだけなので、Rに任せる

 Xs3<-Xs2 %*% t(Xs1)
 Xs4<-rbind(Xs3[1:(N-1),],rep(1,N))
 coefs<-solve(Xs4,c(rep(0,(N-1)),1))
 Apr<-t(Xs1) %*% coefs
 L<-nrm(Apr)
 Apr<-Apr+A
 inside<-FALSE
 if(length(coefs[coefs<0])==0 || (bothside & length(coefs[coefs>0])==0) ){# すべての要素が0以上の条件
  inside<-TRUE
 }
 list(A=A,Xs=Xs,distance=L,Apr=Apr,coefs=coefs,inside=inside)

}

その他、関連ごと（垂線であることのチェック、最短距離であることのチェックなど）

#求めた垂線の足が、AからXs面への最短距離であることを確認
df<-N<-4 # 次元
N<-df-2 # 点の数
#各行が点を表すベクトル
Xs<-matrix(runif(N*df),nrow=N)
#Xs<-matrix(rdirichlet(N,rep(1,df)),nrow=N)
A<-rnorm(df)
L<-MDout2(A,Xs,bothside=FALSE)$distance
testN<-1000
Ls<-rep(0,testN)
tmprs<-matrix(rnorm(testN*N)*100,nrow=testN)
tmpsum<-apply(tmprs,1,sum)
tmprs<-tmprs/tmpsum
for(i in 1:testN){
 tmpA<-t(Xs) %*% tmprs[i,]
 Ls[i]<-sqrt(sum((A-tmpA)^2))
}

plot(sort(Ls),ylim=c(min(L,Ls),max(L,Ls)),type="l")
abline(h=L,col="red")
 min(Ls)-L