Psi | R Documentation |
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.
Psi_mat(data, ind_tri = 0L, use_ind_tri = FALSE, scalar_prod = FALSE,
angles_diff = FALSE)
upper_tri_ind(n)
data |
an array of size |
ind_tri |
if |
use_ind_tri |
use the already computed vector index |
scalar_prod |
return the scalar products
|
angles_diff |
return the (unwrapped) angles difference
|
n |
sample size, used to determine the index vector that gives the
upper triangular part of |
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.
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)
.
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.