How to calculate Pearson correlation coefficient of two 3-dimensional features? The book I read mentions:

Even though, single features do not correlate with each other, sets of features can correlate with other sets of features.

I have an accelerometer, gyroscope, and magnetometer data where all of which are 3-dimensional [x, y, z]. How to calculate if they correlate with each other (for example if accelerometer data correlate with gyroscope data)?
Edit: Adding the exact quote from the book per @whuber request.

However, this pre-processing step may not completely eliminate collinearity, since one or more of the predictors may be functions of two or more of the other predictors.


Second, the univariate importance approach will fail to identify groups of predictors that together have a strong relationship with the response. For example, two predictors may not be highly correlated with the response; however, their interaction may be.


If you are willing to compare the correlation between pairs of vectors, say accelerometer vs gyroscope you could consider computing canonical correlations. The idea is very simple.
Suppose you have blocks $X = (X_1,X_2,X_3)$ (e.g. accelerometer data) and $Y = (Y_1, Y_2, Y_3)$ (gyroscope data), all have the same sample size $n$. Then you consider computing correlations between linear combinations of each, say
$$
r_{ab} = \text{cor}(a^\top X, b^\top Y),
$$
for some $a = (a_1,a_2,a_3)$ and $b = (b_1,b_2,b_3)$. The idea of canonical correlation analysis (CCA) is to look for those $a$ and $b$ which, upon certain constraints, lead to the highest $r_{ab}$ as possible.
If you denote by $S_{1}$ covariance matrix of $X$, $S_2$ the covariance matrix of $Y$ and by $S_{12}$ the cross-covariance of $X,Y$, then the CCA is computed by applying the singular value decomposition (SVD) to the matrix
$$C = S_1^{-1/2} S_{12} S_2^{-1/2}.$$
There can be at most $3$ highest correlations called canonical correlations and these are the singular values of $C$. The first canonical correlation is the most important since it gives the higher correlation possible between the two linear combinations. The optimal vectors $a,b$ are related to the eigenvectors of the SVD of $C$
With CCA you essentially build new (latent) variables by linearly combining the old ones and you try to measure the highest correlation possible between these new variables with the data at hand. If your original variables are not correlated the canonical correlations will be nearly zero. On the other hand, the higher the correlation between the original variables the higher will be the canonical correlations, at least the first canonical correlation.
Canonical correlation analysis can also be seen as a method for dimensionality reduction, e.g. a kind of principal component analysis but, when there are two blocks of variables. If you want to read more about check any textbook on multivariate statistics; e.g. Johnson and Whichern Applied Multivariate Statistical Analysis, Pearson, 6th edition, 2014, ISBN-13: 978-1292024943.
An R implementation of CCA is provided by the CCA package via the cc command or through the cancor command from the stats package.
I'm going to illustrate the use of CCA in a practical example using the dataset iris. In the iris dataset, there are two measures related to sepals, i.e. sepal width and sepal length and two measures related to the petals, petal width and petal length. Suppose we want to have an overall measure of the correlation between sepal measures and petal measures. To do this we run CCA as follows
library(CCA)
X <- iris[,1:2]
Y <- iris[,3:4]
mycca <- cc(X,Y)

> mycca$cor
[1] 0.9409690 0.1239369

The object mycca$cor contains the canonical correlations between the two linear combinations since we have two variables in each block. The first correlation is very high thus we can tell that sepal measures are highly correlated with petal measures.
One might ask if this value is statistically different from zero. To address we can compute a log-likelihood ratio test (again for full details check the aforementioned reference).
rho_hat <- mycca$cor
p = ncol(X)
q = ncol(Y)
n = nrow(X)

# log-likelihood ratio statistic
(w_b = -(n-1-0.5*(p+q+1))*sum(log(1-rho_hat^2)))
[1] 319.6608

# threshold: quantile of the appropriate 
# limiting distribution
> qchisq(p = 0.95, df = p*q)
[1] 9.487729

The test we are conducting here is for $H_0:$ "all canonical correlations are zero" against the alternative $H_1:$ "there is at least one non-zero canonical correlation". Since we reject $H_0$, we conclude that the first canonical correlation is different from zero.
