RiemStiefel (version 0.1.1)

st.pdist: Pairwise Distance for Data on Stiefel Manifold

Description

For data on Stiefel manifold \(x_1,x_2,\ldots,x_N \in St(r,p)\), compute pairwise distances \(d(x_i,x_j)\) via geodesic ("intrinsic") distance or embedded Euclidean "extrinsic" metric. Since computing geodesic has no closed-form expressions, it relies on a numerical approximation which may incur heavier computational burden.

Usage

st.pdist(x, type = c("intrinsic", "extrinsic"), as.dist = TRUE)

Arguments

x

either an array of size \((p\times r\times N)\) or a list of length \(N\) whose elements are \((p\times r)\) matrix on Stiefel manifold.

type

type of distance, either "intrinsic" or "extrinsic".

as.dist

a logical; TRUE to return a dist object or FALSE to return an \((N\times N)\) symmetric matrix.

Value

a dist object or \((N\times N)\) symmetric matrix depending on as.dist.

Examples

# NOT RUN {
#-------------------------------------------------------------------
#      Generate a dataset with two types of Stiefel elements
#-------------------------------------------------------------------
#  group1 : first four columns of (8x8) identity matrix + noise
#  group2 : last  four columns of (8x8) identity matrix + noise

mydata = list()
sdval  = 0.05
diag8  = diag(8)
for (i in 1:10){
  mydata[[i]] = qr.Q(qr(diag8[,1:4] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
}
for (i in 11:20){
  mydata[[i]] = qr.Q(qr(diag8[,5:8] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
}

## compare 'intrinsic' and 'extrinsic' distances
dint = st.pdist(mydata, type="intrinsic", as.dist=FALSE)
dext = st.pdist(mydata, type="extrinsic", as.dist=FALSE)

## visualize
opar = par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
image(dint[,20:1], main="intrinsic")
image(dext[,20:1], main="extrinsic")
par(opar)

# }