Package 'Riemann'

Description

We provide a variety of algorithms for manifold-valued data, including Frechet summaries, hypothesis testing, clustering, visualization, and other learning tasks.

License

Author-owned software in the September 21, 2026 revision of version 0.2.0 is licensed under the GNU General Public License, version 3 (GPL-3). Earlier MIT releases and preserved archives retain their original licenses. Third-party notices are in the installed ‘COPYRIGHTS’ file. Dataset attribution and any data-specific terms are separate; consult the dataset help pages. The separate manuscript and replication workspace has its own license notices.

Author(s)

Maintainer: Kisung You [email protected] (ORCID)

Authors:

• Kisung You [email protected] (ORCID)

See Also

Useful links:

• https://www.kisungyou.com/Riemann/

• Report bugs at https://github.com/kisungyou/Riemann/issues

Description

For a hypersphere $\mathcal{S}^{p-1}$ in $\mathbf{R}^p$, Angular Central Gaussian (ACG) distribution $ACG_p (A)$ is defined via a density

$f(x\vert A) = |A|^{-1/2} (x^\top A^{-1} x)^{-p/2}$

with respect to the uniform measure on $\mathcal{S}^{p-1}$ and $A$ is a symmetric positive-definite matrix. Since $f(x\vert A) = f(-x\vert A)$, it can also be used as an axial distribution on real projective space, which is unit sphere modulo $\lbrace{+1,-1\rbrace}$. One constraint we follow is that $f(x\vert A) = f(x\vert cA)$ for $c > 0$ in that we use a normalized version for numerical stability by restricting $tr(A)=p$.

Usage

dacg(datalist, A, log = FALSE)

racg(n, A)

mle.acg(datalist, ...)

Arguments

Details

The reference measure is the uniform probability measure (total mass one), so the density equals one for the identity parameter. Density evaluation uses Cholesky factors and log determinants; log=TRUE avoids ordinary density overflow. The shape is normalized to trace $p$.

The fixed-point estimate carries convergence, iteration, step, and negative mean-log-likelihood attributes. Failed convergence warns. Existence and uniqueness require sufficient dispersion across all proper linear subspaces; full observed span alone does not establish these statistical conditions. Degenerate updates error, and the implementation does not claim a global existence or uniqueness certificate.

Value

dacg gives a vector of evaluated densities given samples. racg generates unit-norm vectors in $\mathbf{R}^p$ wrapped in a list. mle.acg estimates the SPD matrix $A$.

References

Tyler DE (1987). “Statistical analysis for the angular central Gaussian distribution on the sphere.” Biometrika, 74(3), 579–589. ISSN 0006-3444, 1464-3510.

Mardia KV, Jupp PE (eds.) (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley and Sons, Inc., Hoboken, NJ, USA. ISBN 978-0-470-31697-9 978-0-471-95333-3.

Examples

# -------------------------------------------------------------------
#          Example with Angular Central Gaussian Distribution
#
# Given a fixed A, generate samples and estimate A via ML.
# -------------------------------------------------------------------
## GENERATE AND MLE in R^5
#  Generate data
Atrue = diag(5)          # true SPD matrix
sam1  = racg(50,  Atrue) # random samples
sam2  = racg(100, Atrue)

#  MLE
Amle1 = mle.acg(sam1)
Amle2 = mle.acg(sam2)

#  Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(Atrue[,5:1], axes=FALSE, main="true SPD")
image(Amle1[,5:1], axes=FALSE, main="MLE with n=50")
image(Amle2[,5:1], axes=FALSE, main="MLE with n=100")
par(opar)

Description

As of January 2006, there are 60 cities in the contiguous U.S. with population size larger than $300000$. We extracted information of the cities from the data delivered by maps package. Variables coord and cartesian are two identical representations of locations, which can be mutually converted by sphere.convert.

Usage

data(cities)

Format

a named list containing

names

a length-$60$ vector of city names.

coord

a $(60\times 2)$ matrix of latitude and longitude.

cartesian

a $(60\times 3)$ matrix of cartesian coordinates on the unit sphere.

population

a length-$60$ vector of cities' populations.

See Also

wrap.sphere

Examples

## LOAD THE DATA AND WRAP AS RIEMOBJ
data(cities)
myriem = wrap.sphere(cities$cartesian)

## COMPUTE INTRINSIC/EXTRINSIC MEANS
intmean = as.vector(riem.mean(myriem, geometry="intrinsic")$mean)
extmean = as.vector(riem.mean(myriem, geometry="extrinsic")$mean)

## CONVERT TO GEOGRAPHIC COORDINATES (Lat/Lon)
geo.int = sphere.xyz2geo(intmean)
geo.ext = sphere.xyz2geo(extmean)

Description

Compute density for a fitted mixture model.

Usage

density(object, newdata)

Arguments

Value

a length-$n$ vector of class labels.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

# ---------------------------------------------------- #
#            FIT A MODEL & APPLY THE METHOD
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2    = array(0,c(60,2)) 
for (i in 1:60){
   embed2[i,] = sphere.xyz2geo(locations[i,])
}

# Fit a model
k3 = moSN(locations, k=3)

# Evaluate 
newdensity = density(k3, locations)

Description

A collection of 216 symmetric positive-definite 32-by-32 matrices with four event-label groups. The matrices represent test epochs from one participant in the MNE sample audiovisual experiment, after supervised xDAWN spatial filtering and prototype augmentation. The 32 coordinates comprise 16 class-prototype features and 16 filtered trial features, rather than 32 original EEG or MEG sensors.

Usage

data(ERP)

Format

a named list containing

covariance

an $(32\times 32\times 216)$ array of covariance matrices.

label

a length-$216$ factor with levels LA, LV, RA, RV and counts 57, 57, 53, 49.

Details

A September 21, 2026 reconstruction used the historical pyRiemann 0.2.7 example and numerical routines at commit 176e766f540bd4c7846f38573165fc3d27fc69ca. All 216 stored matrices matched with maximum relative Frobenius error $5.56\times10^{-12}$, and the complete label sequence matched exactly. No rescaling, reordering, or eigenvector-sign correction was needed. The reconstruction used recorded contemporary dependencies; it does not claim byte identity or recovery of the original 2021 runtime.

The recipe reads MNE's prefiltered 0–40 Hz sample recording, sampled at approximately 150.15 Hz after fourfold decimation, applies a 2 Hz IIR high-pass filter, excludes MEG 2443, and selects 305 MEG channels, with EEG excluded. Epochs cover 0–1 second without baseline correction or projection, giving 288 epochs with 151 samples each. A random 25 percent training split (seed 42, without stratification) supplies 72 epochs for supervised xDAWN fitting with four filters per class. The remaining 216 epochs supply the packaged matrices. Covariance estimation uses centered empirical covariance with denominator 151, on the combined prototype and filtered-trial features.

These are repeated epochs from one participant, not independent participants. The representation was learned with labels from the separate 72-epoch training set, so it must not be described as an unsupervised transformation of the original sensors. A later split of the packaged matrices does not recreate raw-data preprocessing or establish performance for new participants.

The matrices retain their original numerical scale; small absolute eigenvalues alone do not imply poor conditioning. No unit conversion, regularization, data values, or labels were changed in this revision. The MNE sample data's public OpenNeuro record, ds000248 version 1.2.4, declares CC0. See the source record for attribution and terms.

Source

Historical pyRiemann ERP recipe, by Pedro Rodrigues (the historical filename says EEG, but the recipe explicitly selects MEG and excludes EEG).

MNE sample data on OpenNeuro, version 1.2.4; pinned source metadata and CC0 declaration. The MNE sample description identifies the audiovisual experiment. The separate writing workspace records the reconstruction in maintenance/reports/erp-provenance-followup-2026-09-21.md.

References

Gramfort A et al. (2013). "MEG and EEG data analysis with MNE-Python." Frontiers in Neuroscience, 7, 267. doi:10.3389/fnins.2013.00267.

Gramfort A et al. (2014). "MNE software for processing MEG and EEG data." NeuroImage, 86, 446–460. doi:10.1016/j.neuroimage.2013.10.027.

See Also

wrap.spd

Examples

## LOAD THE DATA AND WRAP AS RIEMOBJ
data(ERP)
myriem = wrap.spd(ERP$covariance)

Description

This is 29 male and 30 female gorillas' skull landmark data where each individual is represented as 8-ad/landmarks in 2 dimensions. This is a re-arranged version of the data from shapes package.

Usage

data(gorilla)

Format

a named list containing

male

a 3d array of size $(8\times 2\times 29)$

female

a 3d array of size $(8\times 2\times 30)$

Details

These are raw landmark configurations. wrap.landmark removes translation and scale and uses orthogonal alignment, including reflections. The resulting shape analysis is different from analysis of the original sizes or an orientation-preserving quotient. Preserve the documented landmark order in comparisons with shapes.

Source

The arrays match shapes::gorm.dat and shapes::gorf.dat after reordering landmarks by c(1, 5, 4, 3, 2, 8, 7, 6). This exact correspondence was verified against shapes 1.2.8; coordinates and individual ordering are otherwise unchanged. The upstream manual attributes the data to Paul O'Higgins and cites O'Higgins and Dryden (1993). https://cran.r-project.org/package=shapes

References

Dryden IL, Mardia KV (2016). Statistical shape analysis with applications in R, Wiley series in probability and statistics, Second edition edition. John Wiley and Sons, Chichester, UK ; Hoboken, NJ. ISBN 978-1-119-07251-5 978-1-119-07250-8.

O'Higgins P, Dryden IL (1993). "Sexual dimorphism in hominoids: further studies of craniofacial shape differences in Pan, Gorilla, Pongo." Journal of Human Evolution, 24:183–205.

See Also

wrap.landmark

Examples

data(gorilla)                               # load the data
riem.female = wrap.landmark(gorilla$female) # wrap as RIEMOBJ
opar <- par(no.readonly=TRUE)
for (i in 1:30){
  if (i < 2){
    plot(riem.female$data[[i]], cex=0.5, 
         xlim=c(-1,1)/2, ylim=c(-2,2)/5,
         main="30 female gorilla skull preshapes",
         xlab="dimension 1", ylab="dimension 2")
    lines(riem.female$data[[i]])
  } else {
    points(riem.female$data[[i]], cex=0.5)
    lines(riem.female$data[[i]])
  }
}
par(opar)

Description

For a function $f : Gr(k,p) \rightarrow \mathbf{R}$, find the minimizer and the attained minimum value with estimation of distribution algorithm using MACG distribution.

Usage

grassmann.optmacg(func, p, k, ...)

Arguments

Value

a named list containing:

cost

smallest function value attained during the iterative search runs; a global minimum is not guaranteed.

solution

a $(p\times k)$ matrix that attains the cost.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

#-------------------------------------------------------------------
#               Optimization for Eigen-Decomposition
#
# Given (5x5) covariance matrix S, eigendecomposition is can be 
# considered as an optimization on Grassmann manifold. Here, 
# we are trying to find top 3 eigenvalues and compare.
#-------------------------------------------------------------------

## PREPARE
A = cov(matrix(rnorm(100*5), ncol=5)) # define covariance
myfunc <- function(p){                # cost function to minimize
  return(sum(-diag(t(p)%*%A%*%p)))
} 

## SOLVE THE OPTIMIZATION PROBLEM
Aout = grassmann.optmacg(myfunc, p=5, k=3, popsize=100, n.start=30)

## COMPUTE EIGENVALUES
#  1. USE SOLUTIONS TO THE ABOVE OPTIMIZATION 
abase   = Aout$solution
eig3sol = sort(diag(t(abase)%*%A%*%abase), decreasing=TRUE)

#  2. USE BASIC 'EIGEN' FUNCTION
eig3dec = sort(eigen(A)$values, decreasing=TRUE)[1:3]

## VISUALIZE
opar <- par(no.readonly=TRUE)
yran = c(min(min(eig3sol),min(eig3dec))*0.95,
         max(max(eig3sol),max(eig3dec))*1.05)
plot(1:3, eig3sol, type="b", col="red",  pch=19, ylim=yran,
     xlab="index", ylab="eigenvalue", main="compare top 3 eigenvalues")
lines(1:3, eig3dec, type="b", col="blue", pch=19)
legend(1.55, max(yran), legend=c("optimization","decomposition"), col=c("red","blue"),
       lty=rep(1,2), pch=19)
par(opar)

Description

It generates $n$ random samples from Grassmann manifold $Gr(k,p)$.

Usage

grassmann.runif(n, k, p, type = c("list", "array", "riemdata"))

Arguments

Value

an object from one of the above by type option.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

References

Chikuse Y (2003). Statistics on Special Manifolds, volume 174 of Lecture Notes in Statistics. Springer New York, New York, NY. ISBN 978-0-387-00160-9 978-0-387-21540-2.

See Also

stiefel.runif, wrap.grassmann

Examples

#-------------------------------------------------------------------
#                 Draw Samples on Grassmann Manifold 
#-------------------------------------------------------------------
#  Multiple Return Types with 3 Observations of 5-dim subspaces in R^10
dat.list = grassmann.runif(n=3, k=5, p=10, type="list")
dat.arr3 = grassmann.runif(n=3, k=5, p=10, type="array")
dat.riem = grassmann.runif(n=3, k=5, p=10, type="riemdata")

Description

Given the data on Grassmann manifold $Gr(k,p)$, it tests whether the data is distributed uniformly.

Usage

grassmann.utest(grobj, method = c("Bing", "BingM"))

Arguments

Details

This is an asymptotic moment-based test, requiring independent identically distributed observations. It need not detect every nonuniform distribution. The modified expansion is unavailable for ambient dimension two, where its displayed coefficients are singular; full-dimensional subspaces are also excluded. Small-sample calibration and the modified expansion remain experimental; a nominal chi-squared p-value does not certify finite-sample size control.

Value

a (list) object of S3 class htest containing:

statistic

a test statistic.

p.value

$p$-value under $H_0$.

alternative

alternative hypothesis.

method

name of the test.

data.name

name(s) of provided sample data.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

References

Chikuse Y (2003). Statistics on Special Manifolds, volume 174 of Lecture Notes in Statistics. Springer New York, New York, NY. ISBN 978-0-387-00160-9 978-0-387-21540-2.

Mardia KV, Jupp PE (eds.) (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley and Sons, Inc., Hoboken, NJ, USA. ISBN 978-0-470-31697-9 978-0-471-95333-3.

See Also

wrap.grassmann

Examples

#-------------------------------------------------------------------
#   Compare Bingham's original and modified versions of the test
# 
# Test 1. sample uniformly from Gr(2,4)
# Test 2. use perturbed principal components from 'iris' data in R^4
#         which is concentrated around a point to reject H0.
#-------------------------------------------------------------------
## Data Generation
#  1. uniform data
myobj1 = grassmann.runif(n=100, k=2, p=4)

#  2. perturbed principal components
data(iris)
irdat = list()
for (n in 1:100){
   tmpdata    = iris[1:50,1:4] + matrix(rnorm(50*4,sd=0.5),ncol=4)
   irdat[[n]] = eigen(cov(tmpdata))$vectors[,1:2]
}
myobj2 = wrap.grassmann(irdat)

## Test 1 : uniform data
grassmann.utest(myobj1, method="Bing")
grassmann.utest(myobj1, method="BingM")

## Tests : iris data
grassmann.utest(myobj2, method="bINg")   # method names are 
grassmann.utest(myobj2, method="BiNgM")  # CASE - INSENSITIVE !

Description

This dataset contains 10 shapes of 4 subjects's left hands where each shape is represented by 56 landmark points. For each person, first six shapes are equally spaced sequence from maximally to minimally spread fingures. The rest are arbitrarily chosen with two constraints; (1) the palm should face the support and (2) the contour should contain no crossins.

Usage

data(hands)

Format

a named list containing

data

an $(56\times 2\times 40)$ array of landmarks for 40 subjects.

person

a length-$40$ vector of subject indices.

References

Stegmann M, Gomez D (2002) "A Brief Introduction to Statistical Shape Analysis." Informatics and Mathematical Modelling, Technical University of Denmark, DTU.

See Also

wrap.landmark

Examples

## LOAD THE DATA 
data(hands)

## VISUALIZE 6 HANDS OF PERSON 1
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,3))
for (i in 1:6){
  xx = hands$data[,1,i]
  yy = hands$data[,2,i]
  plot(xx,yy,"b", cex=0.9)
}
par(opar)

Description

Given a fitted mixture model of $K$ components, predict labels of observations accordingly.

Usage

label(object, newdata)

Arguments

Value

a length-$n$ vector of class labels.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

# ---------------------------------------------------- #
#            FIT A MODEL & APPLY THE METHOD
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2    = array(0,c(60,2)) 
for (i in 1:60){
   embed2[i,] = sphere.xyz2geo(locations[i,])
}

# Fit a model
k3 = moSN(locations, k=3)

# Evaluate
newlabel = label(k3, locations)

Description

Given a fitted mixture model $f(x)$ and observations $x_1, \ldots, x_n \in \mathcal{M}$, compute the log-likelihood

$L = \log \prod_{i=1}^n f(x_i) = \sum_{i=1}^n \log f(x_i)$

.

Usage

loglkd(object, newdata)

Arguments

Value

the log-likelihood.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

# ---------------------------------------------------- #
#            FIT A MODEL & APPLY THE METHOD
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2    = array(0,c(60,2)) 
for (i in 1:60){
   embed2[i,] = sphere.xyz2geo(locations[i,])
}

# Fit a model
k3 = moSN(locations, k=3)

# Evaluate
newloglkd = round(loglkd(k3, locations), 3)
print(paste0("Log-likelihood for K=3 model fit : ", newloglkd))

Description

Methods for Scalar-Response Kernel Regression Fits

Usage

## S3 method for class 'm2skreg'
fitted(object, ...)

## S3 method for class 'm2skreg'
residuals(object, ...)

## S3 method for class 'm2skreg'
summary(object, ...)

## S3 method for class 'm2skreg'
print(x, ...)

## S3 method for class 'summary.m2skreg'
print(x, ...)

Arguments

Value

fitted and residuals return numeric vectors; residuals are observed responses minus training fitted values. summary returns a summary.m2skreg list. Printing returns its input invisibly.

Description

For Stiefel and Grassmann manifolds $St(r,p)$ and $Gr(r,p)$, the matrix variant of ACG distribution is known as Matrix Angular Central Gaussian (MACG) distribution $MACG_{p,r}(\Sigma)$ with density

$f(X\vert \Sigma) = |\Sigma|^{-r/2} |X^\top \Sigma^{-1} X|^{-p/2}$

where $\Sigma$ is a $(p\times p)$ symmetric positive-definite matrix. Similar to vector-variate ACG case, we follow a convention that $tr(\Sigma)=p$.

Usage

dmacg(datalist, Sigma, log = FALSE)

rmacg(n, r, Sigma)

mle.macg(datalist, ...)

Arguments

Details

The density is relative to the invariant uniform probability measure on the Stiefel manifold, or its push-forward on the Grassmann quotient, and equals one at $\Sigma=I$. Cholesky/log-determinant evaluation avoids determinant products that overflow at extreme parameter scales.

Sampling uses the polar factor of a Gaussian matrix. The fixed-point estimate is trace normalized and carries convergence, iteration, step, and negative mean-log-likelihood attributes. Failed convergence warns; degenerate updates error. Full observed span does not certify the subspace-dispersion conditions needed for an identifiable, unique interior estimate. When $r=p>1$, the density is uniform regardless of $\Sigma$, so estimation errors.

Value

dmacg gives a vector of evaluated densities given samples. rmacg generates $(p\times r)$ orthonormal matrices wrapped in a list. mle.macg estimates the SPD matrix $\Sigma$.

References

Chikuse Y (1990). “The matrix angular central Gaussian distribution.” Journal of Multivariate Analysis, 33(2), 265–274. ISSN 0047259X.

Mardia KV, Jupp PE (eds.) (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley and Sons, Inc., Hoboken, NJ, USA. ISBN 978-0-470-31697-9 978-0-471-95333-3.

Kent JT, Ganeiber AM, Mardia KV (2013). "A new method to simulate the Bingham and related distributions in directional data analysis with applications." arXiv:1310.8110.

See Also

acg

Examples

# -------------------------------------------------------------------
#          Example with Matrix Angular Central Gaussian Distribution
#
# Given a fixed Sigma, generate samples and estimate Sigma via ML.
# -------------------------------------------------------------------
## GENERATE AND MLE in St(2,5)/Gr(2,5)
#  Generate data
Strue = diag(5)                  # true SPD matrix
sam1  = rmacg(n=50,  r=2, Strue) # random samples
sam2  = rmacg(n=100, r=2, Strue) # random samples

#  MLE
Smle1 = mle.macg(sam1)
Smle2 = mle.macg(sam2)

#  Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(Strue[,5:1], axes=FALSE, main="true SPD")
image(Smle1[,5:1], axes=FALSE, main="MLE with n=50")
image(Smle2[,5:1], axes=FALSE, main="MLE with n=100")
par(opar)

Description

For $n$ observations on a $(p-1)$ sphere in $\mathbf{R}^p$, a finite mixture model is fitted whose components are spherical Laplace distributions via the following model

$f(x; \left\lbrace w_k, \mu_k, \sigma_k \right\rbrace_{k=1}^K) = \sum_{k=1}^K w_k SL(x; \mu_k, \sigma_k)$

with parameters $w_k$'s for component weights, $\mu_k$'s for component locations, and $\sigma_k$'s for component scales.

Usage

moSL(
  data,
  k = 2,
  same.sigma = FALSE,
  variants = c("soft", "hard", "stochastic"),
  ...
)

## S3 method for class 'moSL'
loglkd(object, newdata)

## S3 method for class 'moSL'
label(object, newdata)

## S3 method for class 'moSL'
density(object, newdata)

Arguments

Value

a named list of S3 class riemmix containing

cluster

a length-$n$ vector of class labels (from $1:k$).

loglkd

log likelihood of the fitted model.

criteria

a vector of information criteria.

parameters

a list containing proportion, location, and scale. See the section for more details.

membership

an $(n\times k)$ row-stochastic matrix of membership.

Parameters of the fitted model

A fitted model is characterized by three parameters. For $k$-mixture model on a $(p-1)$ sphere in $\mathbf{R}^p$, (1) proportion is a length-$k$ vector of component weight that sums to 1, (2) location is an $(k\times p)$ matrix whose rows are per-cluster locations, and (3) concentration is a length-$k$ vector of scale parameters for each component.

Note on S3 methods

There are three S3 methods; loglkd, label, and density. Given a random sample of size $m$ as newdata, (1) loglkd returns a scalar value of the computed log-likelihood, (2) label returns a length-$m$ vector of cluster assignments, and (3) density evaluates densities of every observation according ot the model fit.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

# ---------------------------------------------------- #
#                 FITTING THE MODEL
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2    = array(0,c(60,2)) 
for (i in 1:60){
   embed2[i,] = sphere.xyz2geo(locations[i,])
}

# Fit the model with different numbers of clusters
k2 = moSL(locations, k=2)
k3 = moSL(locations, k=3)
k4 = moSL(locations, k=4)

# Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3))
plot(embed2, col=k2$cluster, pch=19, main="K=2")
plot(embed2, col=k3$cluster, pch=19, main="K=3")
plot(embed2, col=k4$cluster, pch=19, main="K=4")
par(opar)

# ---------------------------------------------------- #
#                   USE S3 METHODS
# ---------------------------------------------------- #
# Use the same 'locations' data as new data 
# (1) log-likelihood
newloglkd = round(loglkd(k3, locations), 5)
fitloglkd = round(k3$loglkd, 5)
print(paste0("Log-likelihood for K=3 fitted    : ", fitloglkd))
print(paste0("Log-likelihood for K=3 predicted : ", newloglkd))

# (2) label
newlabel = label(k3, locations)

# (3) density
newdensity = density(k3, locations)

Description

For $n$ observations on a $(p-1)$ sphere in $\mathbf{R}^p$, a finite mixture model is fitted whose components are spherical normal distributions via the following model

$f(x; \left\lbrace w_k, \mu_k, \lambda_k \right\rbrace_{k=1}^K) = \sum_{k=1}^K w_k SN(x; \mu_k, \lambda_k)$

with parameters $w_k$'s for component weights, $\mu_k$'s for component locations, and $\lambda_k$'s for component concentrations.

Usage

moSN(
  data,
  k = 2,
  same.lambda = FALSE,
  variants = c("soft", "hard", "stochastic"),
  ...
)

## S3 method for class 'moSN'
loglkd(object, newdata)

## S3 method for class 'moSN'
label(object, newdata)

## S3 method for class 'moSN'
density(object, newdata)

Arguments

Value

a named list of S3 class riemmix containing

cluster

a length-$n$ vector of class labels (from $1:k$).

loglkd

log likelihood of the fitted model.

criteria

a vector of information criteria.

parameters

a list containing proportion, center, and concentration. See the section for more details.

membership

an $(n\times k)$ row-stochastic matrix of membership.

Parameters of the fitted model

A fitted model is characterized by three parameters. For $k$-mixture model on a $(p-1)$ sphere in $\mathbf{R}^p$, (1) proportion is a length-$k$ vector of component weight that sums to 1, (2) center is an $(k\times p)$ matrix whose rows are cluster centers, and (3) concentration is a length-$k$ vector of concentration parameters for each component.

Note on S3 methods

There are three S3 methods; loglkd, label, and density. Given a random sample of size $m$ as newdata, (1) loglkd returns a scalar value of the computed log-likelihood, (2) label returns a length-$m$ vector of cluster assignments, and (3) density evaluates densities of every observation according ot the model fit.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

References

You K, Suh C (2022). “Parameter Estimation and Model-Based Clustering with Spherical Normal Distribution on the Unit Hypersphere.” Computational Statistics and Data Analysis, 107457. ISSN 01679473.

Examples

# ---------------------------------------------------- #
#                 FITTING THE MODEL
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2    = array(0,c(60,2)) 
for (i in 1:60){
   embed2[i,] = sphere.xyz2geo(locations[i,])
}

# Fit the model with different numbers of clusters
k2 = moSN(locations, k=2)
k3 = moSN(locations, k=3)
k4 = moSN(locations, k=4)

# Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3))
plot(embed2, col=k2$cluster, pch=19, main="K=2")
plot(embed2, col=k3$cluster, pch=19, main="K=3")
plot(embed2, col=k4$cluster, pch=19, main="K=4")
par(opar)

# ---------------------------------------------------- #
#                   USE S3 METHODS
# ---------------------------------------------------- #
# Use the same 'locations' data as new data 
# (1) log-likelihood
newloglkd = round(loglkd(k3, locations), 3)
print(paste0("Log-likelihood for K=3 model fit : ", newloglkd))

# (2) label
newlabel = label(k3, locations)

# (3) density
newdensity = density(k3, locations)

Description

The 9 planets in our solar system are evolving the sun via their own orbits. This data provides normal vector of the orbital planes. Normal vectors are unit-norm vectors, so that they are thought to reside on 2-dimensional sphere.

Usage

data(orbital)

Format

an $(9\times 3)$ matrix where each row is a normal vector for a planet.

See Also

wrap.sphere

Examples

## LOAD THE DATA AND WRAP AS RIEMOBJ
data(orbital)
myorb = wrap.sphere(orbital)

## VISUALIZE
mds2d = riem.mds(myorb)$embed
opar <- par(no.readonly=TRUE)
plot(mds2d, main="9 Planets", pch=19, xlab="x", ylab="y")
par(opar)

Description

Passiflora is a genus of about 550 species of flowering plants. This dataset contains 15 landmarks in 2 dimension of 3319 leaves of 40 species. Papers listed in the reference section analyzed the data and found 7 clusters.

Usage

data(passiflora)

Format

a named list containing

data

a 3d array of size $(15\times 2\times 3319)$.

species

a length-$3319$ vector of 40 species factors.

class

a length-$3319$ vector of 7 cluster factors.

References

Chitwood DH, Otoni WC (2017). “Divergent leaf shapes among Passiflora species arise from a shared juvenile morphology.” Plant Direct, 1(5), e00028. ISSN 24754455.

Chitwood DH, Otoni WC (2017). “Morphometric analysis of Passiflora leaves: the relationship between landmarks of the vasculature and elliptical Fourier descriptors of the blade.” GigaScience, 6(1). ISSN 2047-217X.

See Also

wrap.landmark

Examples

data(passiflora)                         # load the data
riemobj = wrap.landmark(passiflora$data) # wrap as RIEMOBJ
pga2d   = riem.pga(riemobj)$embed        # embedding via PGA

opar <- par(no.readonly=TRUE)            # visualize
plot(pga2d, col=passiflora$class, pch=19, cex=0.7,
     main="PGA Embedding of Passiflora Leaves",
     xlab="dimension 1", ylab="dimension 2")
par(opar)

Description

Predicts under the fitted geometry and bandwidth without refitting or using prediction-batch statistics. Training predictors and responses are required.

Usage

## S3 method for class 'm2skreg'
predict(
  object,
  newdata,
  geometry = NULL,
  diagnostics = FALSE,
  block_size = 256L,
  ...
)

Arguments

Value

By default, a numeric vector with one prediction per observation. With diagnostics=TRUE, a list containing prediction and a data frame diagnostics with effective_n, nearest_distance, and nearest_over_bandwidth. The last quantity may be infinite when the ratio is not representable even though the prediction remains finite.

See Also

riem.m2skreg

Description

Standard methods report the observation space, selected geometry, and actual numerical termination. A finite estimate is not by itself proof of convergence. Summary plots display matrix entries or ambient coordinates; they do not imply an isometric view of the manifold.

Usage

## S3 method for class 'riemdata'
print(x, ...)

## S3 method for class 'riemdata'
summary(object, ...)

## S3 method for class 'riem_summary'
print(x, ...)

## S3 method for class 'riem_summary'
summary(object, ...)

## S3 method for class 'riem_summary'
plot(x, ...)

Arguments

Value

Print methods return their object invisibly. Summary methods return structured metadata and numerical diagnostics. Plot methods return invisibly.

Description

Riemann distinguishes validated input representations, available geometry primitives, algorithm contracts, and statistical calibration. A wrapped point or available distance does not establish every algorithm's assumptions.

Details

The installed file ‘method-contracts.csv’ contains one disposition for every exported function and registered S3 method, including its domain, assumptions, restrictions, and evidence scope. Read it with the example below. The four possible dispositions are:

core

The consolidated input, geometry, summary, clustering, tangent PCA, or scalar regression contract, within its documented capabilities.

validated_restricted

A specific legacy contract checked against independent formulas or reference fixtures; the stated restrictions apply.

experimental

A retained legacy interface whose full numerical, mathematical, or statistical contract has not been independently verified.

disabled

An unavailable interface with a documented incompatible implementation. Current geometry restrictions can disable combinations without disabling an entire exported function.

A test-file reference is an evidence pointer, not branch coverage or a proof. No disposition guarantees a global optimum, an identifiable model, finite sample calibration, or validity for every manifold. Experimental methods are retained for compatibility and exploration; their source references and examples should not be read as completed validation.

Common distinctions are particularly relevant: nonnegative graph affinities or regression weights need not form a positive-semidefinite kernel; a metric dissimilarity need not be Euclidean; a local tangent representation need not solve exact nonlinear principal geodesic analysis; and a permutation test requires exchangeability under its specified null. Rayleigh/Bingham and Frechet mean/variance tests need not detect all distributional alternatives.

riem.capabilities() records primitive availability and geometry status. The table here separately records method status. The package maintenance inventory records source and dependency ownership; the release/check logs supply actual execution evidence for each environment.

See Also

riem.capabilities, riem.geometry

Examples

contracts <- read.csv(system.file("method-contracts.csv", package = "Riemann"))
contracts[contracts$name %in% c("riem.kpca", "riem.fanova"),
          c("name", "disposition", "limitations")]

Description

Given $N$ observations $X_1, X_2, \ldots, X_N \in \mathcal{M}$, perform clustering via Competitive Learning Riemannian Quantization (CLRQ). Originally, the algorithm is designed for finding voronoi cells that are used in domain quantization. Given the discrete measure of data, centers of the cells play a role of cluster centers and data are labeled accordingly based on the distance to voronoi centers. For an iterative update of centers, gradient descent algorithm adapted for the Riemannian manifold setting is used with the gain factor sequence

$\gamma_t = \frac{a}{1 + b \sqrt{t}}$

where two parameters $a,b$ are represented by par.a and par.b. For initialization, we provide k-means++ and random seeding options as in k-means.

Usage

riem.clrq(riemobj, k = 2, init = c("plus", "random"), gain.a = 1, gain.b = 1)

Arguments

Value

a named list containing

centers

a 3d array where each slice along 3rd dimension is a matrix representation of class centers.

cluster

a length-$N$ vector of class labels (from $1:k$).

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

References

Le Brigant A, Puechmorel S (2019). “Quantization and clustering on Riemannian manifolds with an application to air traffic analysis.” Journal of Multivariate Analysis, 173, 685–703. ISSN 0047259X.

Bonnabel S (2013). “Stochastic Gradient Descent on Riemannian Manifolds.” IEEE Transactions on Automatic Control, 58(9), 2217–2229. ISSN 0018-9286, 1558-2523.

See Also

riem.kmeans

Examples

#-------------------------------------------------------------------
#          Example on Sphere : a dataset with three types
#
# class 1 : 10 perturbed data points near (1,0,0) on S^2 in R^3
# class 2 : 10 perturbed data points near (0,1,0) on S^2 in R^3
# class 3 : 10 perturbed data points near (0,0,1) on S^2 in R^3
#-------------------------------------------------------------------
## GENERATE DATA
mydata = list()
for (i in 1:10){
  tgt = c(1, stats::rnorm(2, sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 11:20){
  tgt = c(rnorm(1,sd=0.1),1,rnorm(1,sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 21:30){
  tgt = c(stats::rnorm(2, sd=0.1), 1)
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
myriem = wrap.sphere(mydata)
mylabs = rep(c(1,2,3), each=10)

## CLRQ WITH K=2,3,4
clust2 = riem.clrq(myriem, k=2)
clust3 = riem.clrq(myriem, k=3)
clust4 = riem.clrq(myriem, k=4)

## MDS FOR VISUALIZATION
mds2d = riem.mds(myriem, ndim=2)$embed

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2), pty="s")
plot(mds2d, pch=19, main="true label", col=mylabs)
plot(mds2d, pch=19, main="K=2", col=clust2$cluster)
plot(mds2d, pch=19, main="K=3", col=clust3$cluster)
plot(mds2d, pch=19, main="K=4", col=clust4$cluster)
par(opar)

Description

Given manifold-valued data $X_1,X_2,\ldots,X_N \in \mathcal{M}$, this algorithm draws a weighted coreset using the lightweight importance-sampling scheme proposed by the reference below. The Euclidean approximation theorem is not asserted for arbitrary manifold geometries.

Usage

riem.coreset18B(
  riemobj,
  M = max(1L, ceiling(length(riemobj$data)/2)),
  geometry = c("intrinsic", "extrinsic"),
  ...
)

Arguments

Details

Each index is drawn independently with replacement, with probability $q_i = 1/(2N) + d(X_i,\mu)^2/(2\sum_j d(X_j,\mu)^2)$, where $\mu$ is the computed mean in the selected geometry. When all distances are zero the probabilities are uniform. Thus, for any fixed centers, the importance-weighted coreset cost is an unbiased estimate of the full-data squared-distance cost. This identity does not require the computed mean to be globally optimal. Duplicate draws count separately and must not be discarded without summing their weights. Use set.seed() for reproducibility.

Value

a named list containing

coreid

a length-$M$ index vector; repeated indices are retained.

weight

a length-$M$ vector of importance weights, $1/(M q_i)$ for each drawn index $i$.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

References

Bachem O, Lucic M, Krause A (2018). “Scalable k -Means Clustering via Lightweight Coresets.” In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 1119–1127. ISBN 978-1-4503-5552-0.

Examples

#-------------------------------------------------------------------
#          Example on Sphere : a dataset with three types
#
# * 10 perturbed data points near (1,0,0) on S^2 in R^3
# * 10 perturbed data points near (0,1,0) on S^2 in R^3
# * 10 perturbed data points near (0,0,1) on S^2 in R^3
#-------------------------------------------------------------------
## GENERATE DATA
mydata = list()
for (i in 1:10){
  tgt = c(1, stats::rnorm(2, sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 11:20){
  tgt = c(rnorm(1,sd=0.1),1,rnorm(1,sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 21:30){
  tgt = c(stats::rnorm(2, sd=0.1), 1)
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
myriem = wrap.sphere(mydata)

## MDS FOR VISUALIZATION
embed2 = riem.mds(myriem, ndim=2)$embed

## FIND CORESET OF SIZES 3, 6, 9
core1 = riem.coreset18B(myriem, M=3)
core2 = riem.coreset18B(myriem, M=6)
core3 = riem.coreset18B(myriem, M=9)

col1 = rep(1,30); col1[core1$coreid] = 2
col2 = rep(1,30); col2[core2$coreid] = 2
col3 = rep(1,30); col3[core3$coreid] = 2

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
plot(embed2, pch=19, col=col1, main="coreset size=3")
plot(embed2, pch=19, col=col2, main="coreset size=6")
plot(embed2, pch=19, col=col3, main="coreset size=9")
par(opar)

Description

Given two curves $\gamma_1, \gamma_2 : I \rightarrow \mathcal{M}$, we are interested in measuring the discrepancy of two curves. Usually, data are given as discrete observations so we are offering several methods to perform the task. See the section below for detailed description.

Usage

riem.distlp(
  riemobj1,
  riemobj2,
  vect = NULL,
  geometry = c("intrinsic", "extrinsic"),
  ...
)

Arguments

Value

the distance value.

Default Method

Trapezoidal Approximation Assume $\gamma_1 (t_i) = X_i$ and $\gamma_2 (t_i) = Y_i$ for $i=1,2,\ldots,N$. In the Euclidean space, $L_p$ distance between two scalar-valued functions is defined as

$L_p^p (\gamma_1 (x), \gamma_2 (x) = \int_{\mathcal{X}} |\gamma_1 (x) - \gamma_2 (x)|^p dx$

. We extend this approach to manifold-valued curves

$L_p^p (\gamma_1 (t), \gamma_2 (t)) = \int_{t\in I} d^p (\gamma_1 (t), \gamma_2 (t)) dt$

where $d(\cdot,\cdot)$ is an intrinsic/extrinsic distance on manifolds. With the given representations, the above integral is approximated using trapezoidal rule.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

#-------------------------------------------------------------------
#                          Curves on Sphere
#
#  curve1 : y = 0.5*cos(x) on the tangent space at (0,0,1)
#  curve2 : y = 0.5*cos(x) on the tangent space at (0,0,1)
#  curve3 : y = 0.5*sin(x) on the tangent space at (0,0,1)
#
# * distance between curve1 & curve2 should be close to 0.
# * distance between curve1 & curve3 should be large.
#-------------------------------------------------------------------
## GENERATION
vecx  = seq(from=-0.9, to=0.9, length.out=50)
vecy1 = 0.5*cos(vecx) + rnorm(50, sd=0.05)
vecy2 = 0.5*cos(vecx) + rnorm(50, sd=0.05)
vecy3 = 0.5*sin(vecx) + rnorm(50, sd=0.05)

## WRAP AS RIEMOBJ
mat1 = cbind(vecx, vecy1, 1); mat1 = mat1/sqrt(rowSums(mat1^2))
mat2 = cbind(vecx, vecy2, 1); mat2 = mat2/sqrt(rowSums(mat2^2))
mat3 = cbind(vecx, vecy3, 1); mat3 = mat3/sqrt(rowSums(mat3^2))

rcurve1 = wrap.sphere(mat1)
rcurve2 = wrap.sphere(mat2)
rcurve3 = wrap.sphere(mat3)

## COMPUTE DISTANCES
riem.distlp(rcurve1, rcurve2, vect=vecx)
riem.distlp(rcurve1, rcurve3, vect=vecx)

Description

Given two time series - a query $X = (X_1,X_2,\ldots,X_N)$ and a reference $Y = (Y_1,Y_2,\ldots,Y_M)$, riem.dtw computes the most basic version of Dynamic Time Warping (DTW) distance between two series using a symmetric step pattern, meaning no window constraints and others at all. Although the scope of DTW in Euclidean space-valued objects is rich, it is scarce for manifold-valued curves. If you are interested in the topic, we refer to dtw package.

Usage

riem.dtw(riemobj1, riemobj2, geometry = c("intrinsic", "extrinsic"))

Arguments

Value

the distance value.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

#-------------------------------------------------------------------
#                          Curves on Sphere
#
#  curve1 : y = 0.5*cos(x) on the tangent space at (0,0,1)
#  curve2 : y = 0.5*sin(x) on the tangent space at (0,0,1)
# 
#  we will generate two sets for curves of different sizes.
#-------------------------------------------------------------------
## GENERATION
clist = list()
for (i in 1:10){ # curve type 1
  vecx = seq(from=-0.9, to=0.9, length.out=sample(10:50, 1))
  vecy = 0.5*cos(vecx) + rnorm(length(vecx), sd=0.1)
  mats = cbind(vecx, vecy, 1)
  clist[[i]] = wrap.sphere(mats/sqrt(rowSums(mats^2)))
}
for (i in 1:10){ # curve type 2
  vecx = seq(from=-0.9, to=0.9, length.out=sample(10:50, 1))
  vecy = 0.5*sin(vecx) + rnorm(length(vecx), sd=0.1)
  mats = cbind(vecx, vecy, 1)
  clist[[i+10]] = wrap.sphere(mats/sqrt(rowSums(mats^2)))
}

## COMPUTE DISTANCES
outint = array(0,c(20,20))
outext = array(0,c(20,20))
for (i in 1:19){
  for (j in 2:20){
    outint[i,j] <- outint[j,i] <- riem.dtw(clist[[i]], clist[[j]], 
                                           geometry="intrinsic")
    outext[i,j] <- outext[j,i] <- riem.dtw(clist[[i]], clist[[j]],
                                           geometry="extrinsic")
  }
}

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
image(outint[,20:1], axes=FALSE, main="intrinsic DTW Distance")
image(outext[,20:1], axes=FALSE, main="extrinsic DTW Distance")
par(opar)

Description

Compares population Frechet means and variances using the statistic of Dubey and Muller (2019), equation (14). The asymptotic chi-squared calibration targets equality of those summaries, not arbitrary equality of distributions. The permutation version requires exchangeable observations under the stronger null of identical group distributions, and targets alternatives detected by the same mean-and-variance statistic.

Usage

riem.fanova(..., maxiter = 50, eps = 1e-05, geometry = NULL)

riem.fanovaP(..., maxiter = 50, eps = 1e-05, nperm = 99, geometry = NULL)

Arguments

Details

The asymptotic result requires unique population/sample means, positive variances of the squared distances to each group mean, suitable moment/metric-entropy conditions, and group proportions bounded away from zero. The cited paper gives sufficient bounded-space conditions; accepting an input manifold does not establish them for every population distribution. The function rejects failed mean fits and a degenerate empirical denominator.

For group proportions $\lambda_j=n_j/n$, the two denominators are $\sum_j\lambda_j/\widehat\sigma_j^2$ and $\sum_j\lambda_j^2\widehat\sigma_j^2$. Distances are jointly rescaled before forming the statistic; this leaves it unchanged and avoids unit-scale overflow. Numerical rescaling does not transform the underlying geometry.

riem.fanovaP refits all group means for every sampled label assignment. It uses $(1+\#\{T_b\geq T_{obs}\})/(nperm+1)$, including ties. The pooled fit is unchanged by label permutation and is reused. Permutations follow R's current random-number state. The returned Monte Carlo standard error is a plug-in diagnostic for simulation variability, not inferential uncertainty in the scientific effect.

Value

An htest object with statistic, p-value, null/alternative, resolved geometry, calibration, group sizes and mean diagnostics. The permutation version additionally retains permuted statistics, nperm, and mc_se.

References

Dubey P, Müller H (2019). “Fréchet analysis of variance for random objects.” Biometrika, 106(4), 803–821. ISSN 0006-3444, 1464-3510.

Examples

X <- wrap.euclidean(matrix(c(-1, 0, 2, 3), ncol = 1))
Y <- wrap.euclidean(matrix(c(0, 1, 4, 8), ncol = 1))
riem.fanova(X, Y)
set.seed(17)
riem.fanovaP(X, Y, nperm = 19)

Description

Geometry is resolved once when a model is fitted and retained for prediction. For SPD matrices, "intrinsic" means "affine_invariant" and "extrinsic" means "log_euclidean". On the sphere the corresponding names are "round" and "chordal". Euclidean aliases coincide. Landmark geometry currently identifies configurations under the full orthogonal group, including reflections. A capability records an available operation; its validation status must also be considered before interpreting results.

Usage

riem.geometry(riemobj, geometry = NULL)

riem.capabilities()

Arguments

Value

riem.geometry returns a serializable geometry specification. riem.capabilities returns the operation registry as a data frame.

Examples

x <- wrap.spd(list(diag(2), 2 * diag(2)))
riem.geometry(x, "log_euclidean")
riem.capabilities()

Description

Given $N$ observations $X_1, X_2, \ldots, X_M \in \mathcal{M}$, perform hierarchical agglomerative clustering with stats::hclust. The supplied manifold distances are treated as dissimilarities.

Usage

riem.hclust(
  riemobj,
  geometry = NULL,
  method = c("single", "complete", "average", "mcquitty", "ward.D", "ward.D2",
    "centroid", "median"),
  members = NULL
)

Arguments

Details

Ward, centroid and median linkage require care when interpreted as Euclidean sums of squares. In particular, Ward linkage of geodesic distances is not manifold k-means. ward.D is the historical R update and does not implement the Ward criterion implemented by ward.D2. See the stats::hclust documentation for distance powers and member weights.

Value

an object of class hclust. See hclust for details.

References

Müllner D (2013). “fastcluster : Fast Hierarchical, Agglomerative Clustering Routines for R and Python.” Journal of Statistical Software, 53(9). ISSN 1548-7660.

Examples

#-------------------------------------------------------------------
#          Example on Sphere : a dataset with three types
#
# class 1 : 10 perturbed data points near (1,0,0) on S^2 in R^3
# class 2 : 10 perturbed data points near (0,1,0) on S^2 in R^3
# class 3 : 10 perturbed data points near (0,0,1) on S^2 in R^3
#-------------------------------------------------------------------
## GENERATE DATA
mydata = list()
for (i in 1:10){
  tgt = c(1, stats::rnorm(2, sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 11:20){
  tgt = c(rnorm(1,sd=0.1),1,rnorm(1,sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 21:30){
  tgt = c(stats::rnorm(2, sd=0.1), 1)
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
myriem = wrap.sphere(mydata)

## COMPUTE SINGLE AND COMPLETE LINKAGE
hc.sing <- riem.hclust(myriem, method="single")
hc.comp <- riem.hclust(myriem, method="complete")

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
plot(hc.sing, main="single linkage")
plot(hc.comp, main="complete linkage")
par(opar)

Description

Given 2 observations $X_1, X_2 \in \mathcal{M}$, find the interpolated point of a geodesic $\gamma(t)$ for $t \in (0,1)$ which assumes two endpoints $\gamma(0)=X_1$ and $\gamma(1)=X_2$.

Usage

riem.interp(riemobj, t = 0.5, geometry = c("intrinsic", "extrinsic"))

Arguments

Value

an interpolated object in matrix representation on $\mathcal{M}$.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

#-------------------------------------------------------------------
#       Geodesic Interpolation between (1,0) and (0,1) in S^1
#-------------------------------------------------------------------
## PREPARE DATA
sp.start = c(1,0)
sp.end   = c(0,1)
sp.data  = wrap.sphere(rbind(sp.start, sp.end))

## FIND THE INTERPOLATED POINT AT "t=0.25"
mid.int = as.vector(riem.interp(sp.data, t=0.25, geometry="intrinsic"))
mid.ext = as.vector(riem.interp(sp.data, t=0.25, geometry="extrinsic"))

## VISUALIZE
#  Prepare Lines and Points
thetas  = seq(from=0, to=pi/2, length.out=100)
quarter = cbind(cos(thetas), sin(thetas))
pic.pts = rbind(sp.start, mid.int, mid.ext, sp.end)
pic.col = c("black","red","green","black")

# Draw
opar <- par(no.readonly=TRUE)
par(pty="s")
plot(quarter, main="two interpolated points at t=0.25",
     xlab="x", ylab="y", type="l")
points(pic.pts, col=pic.col, pch=19)
text(mid.int[1]-0.1, mid.int[2], "intrinsic", col="red")
text(mid.ext[1]-0.1, mid.ext[2], "extrinsic", col="green")
par(opar)

Description

Given 2 observations $X_1, X_2 \in \mathcal{M}$, find the interpolated points of a geodesic $\gamma(t)$ for $t \in (0,1)$ which assumes two endpoints $\gamma(0)=X_1$ and $\gamma(1)=X_2$.

Usage

riem.interps(
  riemobj,
  vect = c(0.25, 0.5, 0.75),
  geometry = c("intrinsic", "extrinsic")
)

Arguments

Value

a 3d array where $T$ slices along 3rd dimension are interpolated objects in matrix representation.

Validation status

This retained legacy interface is experimental. Its full numerical and statistical contract has not been independently verified across supported inputs. See riem-method-contracts and the installed contract table for method-specific assumptions, restrictions, and evidence scope.

Examples

#-------------------------------------------------------------------
#       Geodesic Interpolation between (1,0) and (0,1) in S^1
#-------------------------------------------------------------------
## PREPARE DATA
sp.start = c(1,0)
sp.end   = c(0,1)
sp.data  = wrap.sphere(rbind(sp.start, sp.end))

## FIND THE INTERPOLATED POINT AT FOR t=0.1, 0.2, ..., 0.9.
myvect  = seq(from=0.1, to=0.9, by=0.1)
geo.int = riem.interps(sp.data, vect=myvect, geometry="intrinsic")
geo.ext = riem.interps(sp.data, vect=myvect, geometry="extrinsic")

geo.int = matrix(geo.int, byrow=TRUE, ncol=2) # re-arrange for plotting
geo.ext = matrix(geo.ext, byrow=TRUE, ncol=2)

## VISUALIZE
#  Prepare Lines and Points
thetas  = seq(from=0, to=pi/2, length.out=100)
quarter = cbind(cos(thetas), sin(thetas))

pts.int = rbind(sp.start, geo.int, sp.end)
pts.ext = rbind(sp.start, geo.ext, sp.end)
col.int = c("black", rep("red",9),  "black")
col.ext = c("black", rep("blue",9), "black")

# Draw
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
plot(quarter, main="intrinsic interpolation", # intrinsic geodesic
     xlab="x", ylab="y", type="l")
points(pts.int, col=col.int, pch=19)
for (i in 1:9){
  text(geo.int[i,1]*0.9, geo.int[i,2]*0.9, 
       paste0(round(i/10,2)), col="red")
}
plot(quarter, main="extrinsic interpolation", # intrinsic geodesic
     xlab="x", ylab="y", type="l")
points(pts.ext, col=col.ext, pch=19)
for (i in 1:9){
  text(geo.ext[i,1]*0.9, geo.ext[i,2]*0.9, 
       paste0(round(i/10,2)), col="blue")
}
par(opar)

Description

ISOMAP - isometric feature mapping - is a dimensionality reduction method to apply classical multidimensional scaling to the geodesic distance that is computed on a weighted nearest neighborhood graph. Nearest neighbor is defined by $k$-NN where two observations are said to be connected when they are mutually included in each other's nearest neighbor. Note that it is possible for geodesic distances to be Inf when nearest neighbor graph construction incurs separate connected components. When an extra parameter padding=TRUE, infinite distances are explicitly replaced by 2 times the maximal finite geodesic distance.

Usage

riem.isomap(riemobj, ndim = 2, nnbd = 5, geometry = NULL, ...)

Arguments

Details

Self-neighbors are excluded and distance ties use observation order. An edge is retained only when each endpoint selects the other among its nearest neighbors. Zero-length edges between duplicate observations are retained. The result records components, adjacency, and whether padding changed disconnected distances. Padding produces an artificial dissimilarity and is not an estimate of the original disconnected graph distance.

Value

a named list containing

embed

an $(N\times ndim)$ matrix whose rows are embedded observations.

References

Silva VD and Tenenbaum JB (2003). "Global Versus Local Methods in Nonlinear Dimensionality Reduction." Advances in Neural Information Processing Systems 15, 721–728. MIT Press.

Examples

#-------------------------------------------------------------------
#          Example on Sphere : a dataset with three types
#
# 10 perturbed data points near (1,0,0) on S^2 in R^3
# 10 perturbed data points near (0,1,0) on S^2 in R^3
# 10 perturbed data points near (0,0,1) on S^2 in R^3
#-------------------------------------------------------------------
## GENERATE DATA
mydata = list()
for (i in 1:10){
  tgt = c(1, stats::rnorm(2, sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 11:20){
  tgt = c(rnorm(1,sd=0.1),1,rnorm(1,sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 21:30){
  tgt = c(stats::rnorm(2, sd=0.1), 1)
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
myriem = wrap.sphere(mydata)
mylabs = rep(c(1,2,3), each=10)

## MDS AND ISOMAP WITH DIFFERENT NEIGHBORHOOD SIZE
mdss = riem.mds(myriem)$embed
iso1 = riem.isomap(myriem, nnbd=5, padding=TRUE)$embed
iso2 = riem.isomap(myriem, nnbd=10, padding=TRUE)$embed

## VISUALIZE
opar = par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
plot(mdss, col=mylabs, pch=19, main="MDS")
plot(iso1, col=mylabs, pch=19, main="ISOMAP:nnbd=5")
plot(iso2, col=mylabs, pch=19, main="ISOMAP:nnbd=10")
par(opar)

Description

Given $N$ observations $X_1, X_2, \ldots, X_N \in \mathcal{M}$, perform k-means clustering by minimizing within-cluster sum of squares (WCSS). Since the problem is NP-hard and sensitive to the initialization, we provide an option with multiple starts and return the best result with respect to WCSS.

Usage

riem.kmeans(riemobj, k = 2, geometry = NULL, ...)

## S3 method for class 'riem_kmeans'
predict(object, newdata, type = c("class", "distance"), ...)

## S3 method for class 'riem_kmeans'
print(x, ...)

## S3 method for class 'riem_kmeans'
summary(object, ...)

Arguments

Details

The default is now Lloyd iteration, with each cluster updated by the shared Fréchet-mean solver. A Lloyd objective increase beyond tolerance is rejected and reported. The retained MacQueen option sequentially assigns observations and recomputes affected cluster means; no global monotonicity guarantee is claimed for approximate or nonconvex manifold means. Tied distances choose the lowest cluster index. Empty clusters are reseeded using a maximal-residual observation from a nonsingleton cluster; impossible separation and failed mean solves invalidate that start. Returned labels always assign observations to the final stored centers. Nonconverged finite starts remain visible, and returning one produces a warning. Use set.seed() to reproduce random initializations.

Value

a fitted riem_kmeans object containing

cluster

a length-$N$ vector of class labels (from $1:k$).

means

a 3d array where each slice along 3rd dimension is a matrix representation of class mean.

score

within-cluster sum of squares (WCSS).

References

Lloyd S (1982). “Least squares quantization in PCM.” IEEE Transactions on Information Theory, 28(2), 129–137. ISSN 0018-9448.

MacQueen J (1967). “Some methods for classification and analysis of multivariate observations.” In Proceedings of the fifth berkeley symposium on mathematical statistics and probability, volume 1: Statistics, 281–297.

See Also

riem.kmeanspp

Examples

#-------------------------------------------------------------------
#          Example on Sphere : a dataset with three types
#
# class 1 : 10 perturbed data points near (1,0,0) on S^2 in R^3
# class 2 : 10 perturbed data points near (0,1,0) on S^2 in R^3
# class 3 : 10 perturbed data points near (0,0,1) on S^2 in R^3
#-------------------------------------------------------------------
## GENERATE DATA
mydata = list()
for (i in 1:10){
  tgt = c(1, stats::rnorm(2, sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 11:20){
  tgt = c(rnorm(1,sd=0.1),1,rnorm(1,sd=0.1))
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
for (i in 21:30){
  tgt = c(stats::rnorm(2, sd=0.1), 1)
  mydata[[i]] = tgt/sqrt(sum(tgt^2))
}
myriem = wrap.sphere(mydata)
mylabs = rep(c(1,2,3), each=10)

## K-MEANS WITH K=2,3,4
clust2 = riem.kmeans(myriem, k=2)
clust3 = riem.kmeans(myriem, k=3)
clust4 = riem.kmeans(myriem, k=4)

## MDS FOR VISUALIZATION
mds2d = riem.mds(myriem, ndim=2)$embed

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2), pty="s")
plot(mds2d, pch=19, main="true label", col=mylabs)
plot(mds2d, pch=19, main="K=2", col=clust2$cluster)
plot(mds2d, pch=19, main="K=3", col=clust3$cluster)
plot(mds2d, pch=19, main="K=4", col=clust4$cluster)
par(opar)

Description

Apply weighted Lloyd iterations to an independently sampled lightweight coreset of manifold-valued observations. The Euclidean coreset approximation theorem is not asserted for arbitrary manifold geometries.

Usage

riem.kmeans18B(
  riemobj,
  k = 2,
  M = max(k, ceiling(length(riemobj$data)/2)),
  geometry = c("intrinsic", "extrinsic"),
  ...
)

Arguments