2017-07-27

6　Kendall's (Direct) Similarity Shape Spaces ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

d次元球面上の点集合があったときに、m=2,3次元回転を同一視して、商群を取ったもの
$S^d/G$
$G = SO(m)$ : 特殊直交群 (mxm直交行列であって、行列式が +1)
m=2,3が大事な例
特徴点がm=2次元座標として得られているとする
特徴点の数k=30の例
まず、m=2次元の両軸について原点中心にする(自由度が2減る)
さらに、mxk次元ベクトルとして、「長さ」を１に標準化する
自由度はmkから(2+1)下がってmk-3

m <- 2
k <- 30

z <- matrix(rnorm(m*k),nrow=m)
z[1,] <- z[1,]*0.3 + 1
z[2,] <- z[2,]*0.2 + 2

ang <- Arg(z[1,] + i * z[2,])
z <- z[,order(ang)]


plot(t(z),pch=20,asp=1)

z.normed <- z - apply(z,1,mean)
plot(t(z.normed),pch=20,asp=1)

z.normed. <- z.normed/sqrt(sum(z.normed^2))

plot(t(cbind(z,z.normed,z.normed.)),pch=20,asp=1,type="p")
points(t(z),pch=20,col=1,type="b")
points(t(z.normed),pch=20,col=2,type="b")
points(t(z.normed.),pch=20,col=3,type="b"

自由度mxkをmxk- m - 1に下げた。そのうち、最後の「ノルム=1」の自由度の減は、線型変換ではないので、以下に示すように、線型変換によってmxk-m-1次元平面に映すことは出来ないが、mxk-m次元平面にある、単位球に移すことは出来る
- 以下では、多数のm=2,k=3の長さmxk=6のベクトルを作成し、それを重心標準化、ノルム標準化することによって得られる、6次元空間座標を持った多様体上の点が、mxk-m=4次元平坦空間に存在し、それがノルム=1の制約を持っていること、すなわち、mxk-m次元にある $S^{mxk-m-1}$ という多様体を成していることを示すRのコードである

library(MASS)
library(Matrix)
library(rgl)

n <- 1000 # データセットの数
# そのmxn次元座標を納める行列
Z <- matrix(0,n,m*k)

for(i in 1:n){
	z <- matrix(rnorm(m*k),nrow=m)
	z[1,] <- z[1,]*0.3 + 1
	z[2,] <- z[2,]*0.2 + 2
	z.normed <- z - apply(z,1,mean)
	z.normed. <- z.normed/sqrt(sum(z.normed^2))
	Z[i,] <- z.normed.
}

# ランクがmxk-m
# mxk-m次元空間に収まっていることがわかる
rankMatrix(Z)
# QR分解を用いて、mxk-m次元空間に移す回転行列を計算する
# 発生させたmxkベクトルのうち、ランク本のベクトルを用いる
df <- m*k - (m)
topZ <- Z[1:df,]
# QR分解
qr.out <- qr(t(topZ))
Q <- qr.Q(qr.out)
R <- qr.R(qr.out)

# QR分解により、t(topZ)のペアワイズなベクトル内積が
# R(減次元ベクトルの束)のペアワイズなベクトル内積と一致することを確認する
# 以下はessentially zero
range(t(R) %*% R - topZ %*% t(topZ))

# t(topZ) = QR と出来たから
# Q.inv t(topZ) = R のQ.invを作る
Q.inv <- ginv(Q)
# 検算
round(Q.inv %*% t(topZ) - R)
# Q.invを使って、すべてのmxkベクトルを減次元写像する
rotZ <- Q.inv %*% t(Z)
# ペアワイズ内積が保たれていることを確認する
d1 <- Z %*% t(Z)
d2 <- t(rotZ) %*% rotZ
range(d1-d2)

# 減次元座標のノルムの確認
range(apply(rotZ^2,2,sum))

# mxk -m = 4次元空間の球面はなので、３次元空間にプロットすると
# 球面と球の内部に点が分布する
plot3d(t(rotZ[1:3,]))
plot3d(t(rotZ[2:4,]))

> rankMatrix(Z)
[1] 4
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 2.220446e-13
> # QR分解により、t(topZ)のペアワイズなベクトル内積が
> # R(減次元ベクトルの束)のペアワイズなベクトル内積と一致することを確認する
> # 以下はessentially zero
> range(t(R) %*% R - topZ %*% t(topZ))
[1] -1.110223e-16  1.110223e-16
> # 検算
> round(Q.inv %*% t(topZ) - R)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0
> range(d1-d2)
[1] -1.221245e-15  1.221245e-15
> 
> # 減次元座標のノルムの確認
> range(apply(rotZ^2,2,sum))
[1] 1 1

この $S^{m\times k - m -1}$ 球面は、Preshape sphereと呼ぶ。なぜなら、まだ回転同一視をしていないから
m=2の場合、回転は群 SO(2)であってそれは $S^1$ (円環)なので、 $S^{m\times k - m- 1}/G = S^{2k-3}/S^1$ なる、2k-4次元のコンパクトな多様体であることがわかる
これはk-2次元のcomplex projective spaceと(微分)同相
これは次のようにも書ける
- - 行列z、ただし、Traceの制約は、zの全要素の二乗の和が1であること、1_kなる、全成分が1の長さkのベクトルを含む式は、k要素のすべての平均が0でるように標準化した制約に相当する
- うまく等長変換して、mxk-m次元平面に移せるので
  - $S_m^k = \{z \in R^{m \times k -m} | Trace(z t(z)) =1\}$ とまで持ち込める
  - こうすることで、線型制約が無くなるので便利
  - 本ではHelmert変換を使って-m次元する、と書いている
  - 上のコードでは、QR分解を使っている
さて、回転同一視をする
Direct Similarity Shape Spaces
- $\Sigma_m^k = S_m^k/SO(m)$
- ではは多様体になるとは限らない
  - $m\times k$ 行列のランクがm-1以上なら、多様体になる
- そのような場合を $NS_m^k = \{z \in S_m^k|rank(z) \ge m-1\}$ と書けば
- $\Sigma_{0m}^k = NS_m^k/SO(m)$ で、これは多様体であって、その次元は $km-m-1-\frac{m(m-1)}{2}$ だという
- 微分同相が担保されれば、接空間、接空間を軌道とその直交成分に分けることなどの議論が進められる
- ※　行列と転置行列の積のトレースは、行列をベクトル化して、その内積を取るのと同じであることに注意
- 接空間は $T_z S_m^k = \{v \in H_m^k| Trace (v t(z))=0\}$ 、ただし $H_m^k$ は $S_m^k$ を含んだ超平面。その超平面内の球上の接空間は球面上のベクトル $z$ との直交条件で示せる
回転同一視をして、その上での形間(測地線)距離は $min_{T \in SO(m)} d_{gs}(x,Ty)$ 、ただし $d_{gs}(x,y) = arccos(Trace(yt(x))$ と表される、球面上の測地線
したがって $arccos(max_{T \in SO(m)} Trace(Ty, t(x)))$
の最大化はのsvd分解によって求まる
- 最大値 $y x^T = U\Lambda V$ というsvd分解での $Trace(\Lambda) = \sum \lambda_i$ であり
- それをもたらす回転は $T= U^T V^T$

library(svd)

# x,y,z座標球面スカラー場に対して
# k個のスペクトルが得られたとする
# それが形情報
# x,y なるm * k 行列

# 3次元ベクトルがk個
m <- 3
k <- 5

# x,y,z座標を平均0にして
# 全要素の二乗和を1に標準化する
x <- matrix(rnorm(m*k),ncol=k)
y <- matrix(rnorm(m*k),ncol=k)

x <- t(t(x)-apply(x,2,mean))
y <- t(t(y)-apply(y,2,mean))
x <- x/sqrt(sum(x^2))
y <- y/sqrt(sum(y^2))

# 検算
apply(x,1,sum)
apply(y,1,sum)
sum(x^2)
sum(y^2)

# svd分解法
# http://www2.stat.duke.edu/~ab216/monograph.pdf の
# p84 の式(6.4)の辺りを使って
# yを回転して、xと近づける最適化を行う

yxt <- y %*% t(x)

svdout <- svd(yxt)
U <- svdout$u
V <- svdout$v

# svd 分解の確認
# テキストの式表現とsvd()関数の式表現のずれに注意する
U %*% diag(svdout$d) %*% t(V)
yxt

# テキストで言うところの最適回転行列 TをWとする
W <- t(t(V)) %*% t(U)

# Trace(Ty t(x)) = sum(eigenvalue)の検算

sum(diag(W %*% y %*% t(x)))
sum(svdout$d)
W %*% y