# R/gr_pdist2.R In RiemGrassmann: Inference, Learning, and Optimization on Grassmann Manifold

#### Documented in gr.pdist2

#' Pairwise Distance for Two Sets Data on Grassmann Manifold
#'
#' For data on grassmann manifold \eqn{x_1,x_2,\ldots,x_M \in Gr(k,n)} and
#' \eqn{y_1,y_2,\ldots,y_N \in Gr(k,n)}, compute pairwise distances \eqn{d(x_i,y_j)} via several metrics. The distance type \code{"intrinsic"} corresponds to
#' geodesic distance while \code{"extrinsic"} is equivalent to \code{"chordal"} distance.
#'
#'
#' @param input1 either an array of size \eqn{(n\times k\times M)} or a list of length \eqn{M} whose elements are \eqn{(n\times k)} orthonormal basis (ONB) on Grassmann manifold.
#' @param input2 either an array of size \eqn{(n\times k\times N)} or a list of length \eqn{N} whose elements are \eqn{(n\times k)} orthonormal basis (ONB) on Grassmann manifold.
#' @param type type of distance measure. Name of each type is \emph{Case Insensitive} and \emph{hyphen} can be omitted.
#' @param useR a logical; \code{TRUE} to use R computations while \code{FALSE} goes everything in C++.
#'
#' @return an \eqn{(M\times N)} matrix of pairwise distances.
#'
#' @examples
#' ## generate a dataset with two types of Grassmann elements
#' #  group1 : first four columns of (8x8) identity matrix + noise
#' #  group2 : last  four columns of (8x8) identity matrix + noise
#'
#' mydata1 = list()
#' mydata2 = list()
#' sdval   = 0.25
#' diag8   = diag(8)
#' for (i in 1:10){
#'   mydata1[[i]] = qr.Q(qr(diag8[,1:4] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
#' }
#' for (i in 1:10){
#'   mydata2[[i]] = qr.Q(qr(diag8[,5:8] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
#' }
#'
#'
#' \donttest{
#' ## compute and visualize distances for all types
#' #  we will iterate over all measures
#' alltypes = c("intrinsic","extrinsic","asimov","binetcauchy",
#' "chordal","fubinistudy","martin","procrustes","projection","spectral")
#' ntypes   = length(alltypes)
#'
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(2,5), pty="s")
#' for (i in 1:ntypes){
#'   dmat = gr.pdist2(mydata1, mydata2, type=alltypes[i])
#'   image(dmat, main=alltypes[i])
#' }
#' par(opar)
#' }
#'
#' @author Kisung You
#' @export
gr.pdist2 <- function(input1, input2, type=c("Intrinsic","Extrinsic","Asimov","Binet-Cauchy",
"Chordal","Fubini-Study","Martin","Procrustes",
"Projection","Spectral"), useR=FALSE){
############################################################
# Preprocessing
input1 = RiemBase::riemfactory(return_gr(input1), name="grassmann")
input2 = RiemBase::riemfactory(return_gr(input2), name="grassmann")
alltypes = c("intrinsic","extrinsic","asimov","binetcauchy",
"chordal","fubinistudy","martin","procrustes",
"projection","spectral")
if (missing(type)){
mytype = "intrinsic"
} else {
mytype = match.arg(gsub("[-]","",tolower(type)),alltypes)
}
myuseR = as.logical(useR)

############################################################
# Computation
if ((all(mytype=="intrinsic"))||(all(mytype=="extrinsic"))){
output = RiemBaseExt::rstat.pdist2(input1, input2, type=mytype)
} else {
if (myuseR){
output = gr.pdist2.others(input1, input2, method=mytype)
} else {
output = cpp_pdist2(input1$data, input2$data, mytype)
}
}

############################################################
# Return
return(output)
}

# auxiliary here ----------------------------------------------------------
#' @keywords internal
#' @noRd
gr.pdist2.others <- function(input1, input2, method){
gdat1 = input1$data gdat2 = input2$data

N = length(gdat1)
M = length(gdat2)

output = array(0,c(N,M))
for (n in 1:N){
Xn = as.matrix(gdat1[[n]])
for (m in 1:M){
Xm = as.matrix(gdat2[[m]])
output[n,m] <- cdist_distance(Xn,Xm,method)
}
}
return(output)
}

# ## personal test
# mydata1 = list()
# mydata2 = list()
# sdval   = 0.25
# diag8   = diag(8)
# for (i in 1:10){
#   mydata1[[i]] = qr.Q(qr(diag8[,1:4] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
# }
# for (i in 1:10){
#   mydata2[[i]] = qr.Q(qr(diag8[,5:8] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
# }
#
# alltypes = c("intrinsic","extrinsic","asimov","binet-cauchy",
# "chordal","fubini-study","martin","procrustes","projection","spectral")
# ntypes   = length(alltypes)
# for (i in 1:ntypes){
#   dist.R = gr.pdist2(mydata1, mydata2, type=alltypes[i], useR=TRUE)
#   dist.C = gr.pdist2(mydata1, mydata2, type=alltypes[i], useR=FALSE)
#   print(paste0(alltypes[i]," difference=",norm(dist.R-dist.C,"F")))
# }
#
# library(microbenchmark)
# microbenchmark::microbenchmark(
#   "R" = gr.pdist2(mydata1, mydata2, type=alltypes[i], useR=TRUE),
#   "C" = gr.pdist2(mydata1, mydata2, type=alltypes[i], useR=FALSE)
# )


## Try the RiemGrassmann package in your browser

Any scripts or data that you put into this service are public.

RiemGrassmann documentation built on March 25, 2020, 5:07 p.m.