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

#### Documented in gr.pdist

#' Pairwise Distance for Data on Grassmann Manifold
#'
#' For data on grassmann manifold \eqn{x_1,x_2,\ldots,x_N \in Gr(k,n)}, compute pairwise distances \eqn{d(x_i,x_j)} via several metrics. The distance type \code{"intrinsic"} corresponds to
#' geodesic distance while \code{"extrinsic"} is equivalent to \code{"chordal"} distance.
#'
#' @param input 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 as.dist a logical; \code{TRUE} to return a \code{\link[stats]{dist}} object or \code{FALSE} to return an \eqn{(N\times N)} symmetric matrix.
#' @param useR a logical; \code{TRUE} to use R computations while \code{FALSE} goes everything in C++.
#'
#' @return a \code{\link[stats]{dist}} object or \eqn{(N\times N)} symmetric matrix depending on \code{as.dist}.
#'
#' @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
#'
#' mydata = list()
#' sdval  = 0.25
#' 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)))
#' }
#'
#' ## try 'intrinsic' distance using C++ implementation
#' dmat = gr.pdist(mydata, type="intrinsic", as.dist=FALSE)
#' opar = par(no.readonly=TRUE)
#' par(pty="s")
#' image(dmat, main="intrinsic distance")
#' par(opar)
#'
#' \donttest{
#' ## compute and visualize distances for all types
#' #  we will iterate over all measures
#' alltypes = c("intrinsic","extrinsic","asimov","binet-cauchy",
#' "chordal","fubini-study","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.pdist(mydata, type=alltypes[i], as.dist=FALSE)
#'   image(dmat[,20:1], main=alltypes[i])
#' }
#' par(opar)
#' }
#'
#' @author Kisung You
#' @export
gr.pdist <- function(input, type=c("Intrinsic","Extrinsic","Asimov","Binet-Cauchy",
"Chordal","Fubini-Study","Martin","Procrustes",
"Projection","Spectral"), as.dist=TRUE, useR=FALSE){
############################################################
# Preprocessing
x        = RiemBase::riemfactory(return_gr(input), 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)
}
retdist = as.logical(as.dist)
myuseR  = as.logical(useR)

############################################################
# Computation
if ((all(mytype=="intrinsic"))||(all(mytype=="extrinsic"))){
output = RiemBaseExt::rstat.pdist(x, type=mytype, as.dist=FALSE)
} else {
if (myuseR){
output = gr.pdist.others(x, method=mytype)
} else {
output = cpp_pdist(x$data, mytype) } } ############################################################ # Return if (retdist){ return(stats::as.dist(output)) } else { return(output) } } # auxiliary here ---------------------------------------------------------- #' @keywords internal #' @noRd gr.pdist.others <- function(input, method){ grdat = input$data
N     = length(grdat)

output = array(0,c(N,N))
for (i in 1:(N-1)){
Xi = as.matrix(grdat[[i]])
for (j in (i+1):N){
Xj = as.matrix(grdat[[j]])
output[i,j] <- output[j,i] <- cdist_distance(Xi,Xj,method)
}
}
return(output)
}
#' @keywords internal
#' @noRd
gr.pdist.nocheck <- function(input, mytype, asdist=TRUE){
if ((all(mytype=="intrinsic"))||(all(mytype=="extrinsic"))){
output = RiemBaseExt::rstat.pdist(input, type=mytype, as.dist=FALSE)
} else {
output = gr.pdist.others(input, method=mytype)
}
if (asdist){
return(stats::as.dist(output))
} else {
return(output)
}
}

# ## difference between C and R
# mydata = list()
# sdval  = 0.25
# 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)))
# }
#
# #  we will iterate over all measures
# 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.pdist(mydata, type=alltypes[i], useR=TRUE, as.dist=FALSE)
#   dist.C = gr.pdist(mydata, type=alltypes[i], useR=FALSE, as.dist=FALSE)
#   print(paste0(alltypes[i]," difference=",norm(dist.R-dist.C,"F")))
# }
#
# library(microbenchmark)
# microbenchmark::microbenchmark(
#   "R" = gr.pdist(mydata, type=alltypes[i], useR=TRUE, as.dist=FALSE),
#   "C" = gr.pdist(mydata, type=alltypes[i], useR=FALSE, as.dist=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.