Psi {sphunif}

Shortest angles matrix

Description

Efficient computation of the shortest angles matrix \boldsymbol\Psi, defined as

\Psi_{ij}:=\cos^{-1}({\bf X}_i'{\bf X}_j),\quad i,j=1,\ldots,n,

for a sample {\bf X}_1,\ldots,{\bf X}_n\in S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}, p\ge 2.

For a circular sample \Theta_1, \ldots, \Theta_n \in [0, 2\pi), \boldsymbol\Psi can be expressed as

\Psi_{ij}=\pi-|\pi-|\Theta_i-\Theta_j||,\quad i,j=1,\ldots,n.

Usage

Psi_mat(data, ind_tri = 0L, use_ind_tri = FALSE, scalar_prod = FALSE,
  angles_diff = FALSE)

upper_tri_ind(n)

Arguments

data

an array of size c(n, p, M) containing the Cartesian coordinates of M samples of size n of directions on S^{p-1}. Alternatively if p = 2, an array of size c(n, 1, M) containing the angles on [0, 2\pi) of the M circular samples of size n on S^{1}. Must not contain NA's.

ind_tri

if use_ind_tri = TRUE, the vector of 0-based indexes provided by upper_tri_ind(n), which allows to extract the upper triangular part of the matrix \boldsymbol\Psi. See the examples.

use_ind_tri

use the already computed vector index ind_tri? If FALSE (default), ind_tri is computed internally.

scalar_prod

return the scalar products {\bf X}_i'{\bf X} instead of the shortest angles? Only taken into account for data in Cartesian form. Defaults to FALSE.

angles_diff

return the (unwrapped) angles difference \Theta_i-\Theta_j instead of the shortest angles? Only taken into account for data in angular form. Defaults to FALSE.

n

sample size, used to determine the index vector that gives the upper triangular part of \boldsymbol\Psi.

Value

• Psi_mat: a matrix of size c(n * (n - 1) / 2, M) containing, for each column, the vector half of \boldsymbol\Psi for each of the M samples.

• upper_tri_ind: a matrix of size n * (n - 1) / 2 containing the 0-based linear indexes for extracting the upper triangular matrix of a matrix of size c(n, n), diagonal excluded, assuming column-major order.

Warning

Be careful on avoiding the next bad usages of Psi_mat, which will produce spurious results:

• The directions in data do not have unit norm when Cartesian coordinates are employed.

• The entries of data are not in [0, 2\pi) when polar coordinates are employed.

• ind_tri is a vector of size n * (n - 1) / 2 that does not contain the indexes produced by upper_tri_ind(n).

Examples

# Shortest angles
n <- 5
X <- r_unif_sph(n = n, p = 2, M = 2)
Theta <- X_to_Theta(X)
dim(Theta) <- c(n, 1, 2)
Psi_mat(X)
Psi_mat(Theta)

# Precompute ind_tri
ind_tri <- upper_tri_ind(n)
Psi_mat(X, ind_tri = ind_tri, use_ind_tri = TRUE)

# Compare with R
A <- acos(tcrossprod(X[, , 1]))
ind <- upper.tri(A)
A[ind]

# Reconstruct matrix
Psi_vec <- Psi_mat(Theta[, , 1, drop = FALSE])
Psi <- matrix(0, nrow = n, ncol = n)
Psi[upper.tri(Psi)] <- Psi_vec
Psi <- Psi + t(Psi)