For each column \(j=1,\ldots,p\), let
\(R_{\cdot j} \in \{1,\ldots,n\}^n\) denote the (mid-)ranks of
\(X_{\cdot j}\), assigning average ranks to ties. The mean rank is
\(\bar R_j = (n+1)/2\) regardless of ties. Define the centred rank
vectors \(\tilde R_{\cdot j} = R_{\cdot j} - \bar R_j \mathbf{1}\),
where \(\mathbf{1}\in\mathbb{R}^n\) is the all-ones vector. The
Spearman correlation between columns \(i\) and \(j\) is the Pearson
correlation of their rank vectors:
$$
\rho_S(i,j) \;=\;
\frac{\sum_{k=1}^n (R_{ki}-\bar R_i)(R_{kj}-\bar R_j)}
{\sqrt{\sum_{k=1}^n (R_{ki}-\bar R_i)^2}\;
\sqrt{\sum_{k=1}^n (R_{kj}-\bar R_j)^2}}.
$$
In matrix form, with \(R=[R_{\cdot 1},\ldots,R_{\cdot p}]\),
\(\mu=(n+1)\mathbf{1}_p/2\) for \(\mathbf{1}_p\in\mathbb{R}^p\), and
\(S_R=\bigl(R-\mathbf{1}\mu^\top\bigr)^\top
\bigl(R-\mathbf{1}\mu^\top\bigr)/(n-1)\),
the Spearman correlation matrix is
$$
\widehat{\rho}_S \;=\; D^{-1/2} S_R D^{-1/2}, \qquad
D \;=\; \mathrm{diag}(\mathrm{diag}(S_R)).
$$
When there are no ties, the familiar rank-difference formula obtains
$$
\rho_S(i,j) \;=\; 1 - \frac{6}{n(n^2-1)} \sum_{k=1}^n d_k^2,
\quad d_k \;=\; R_{ki}-R_{kj},
$$
but this expression does not hold under ties; computing Pearson on
mid-ranks (as above) is the standard tie-robust approach. Without ties,
\(\mathrm{Var}(R_{\cdot j})=(n^2-1)/12\); with ties, the variance is
smaller.
\(\rho_S(i,j) \in [-1,1]\) and \(\widehat{\rho}_S\) is symmetric
positive semi-definite by construction (up to floating-point error). The
implementation symmetrises the result to remove round-off asymmetry.
Spearman's correlation is invariant to strictly monotone transformations
applied separately to each variable.
Computation. Each column is ranked (mid-ranks) to form \(R\).
The product \(R^\top R\) is computed via a 'BLAS' symmetric rank update
('SYRK'), and centred using
$$
(R-\mathbf{1}\mu^\top)^\top (R-\mathbf{1}\mu^\top)
\;=\; R^\top R \;-\; n\,\mu\mu^\top,
$$
avoiding an explicit centred copy. Division by \(n-1\) yields the sample
covariance of ranks; standardising by \(D^{-1/2}\) gives \(\widehat{\rho}_S\).
Columns with zero rank variance (all values equal) are returned as NA
along their row/column; the corresponding diagonal entry is also NA.
When check_na = FALSE, each \((i,j)\) estimate is recomputed on the
pairwise complete-case overlap of columns \(i\) and \(j\). When
ci = TRUE, confidence intervals are computed in 'C++' using the
jackknife Euclidean-likelihood method of de Carvalho and Marques (2012).
For a pairwise estimate \(U = \hat\rho_S\), delete-one jackknife
pseudo-values are formed as
$$
Z_i = nU - (n-1)U_{(-i)}, \qquad i = 1,\ldots,n,
$$
where \(U_{(-i)}\) is the Spearman correlation after removing observation
\(i\). The confidence limits solve
$$
\frac{n(U-\theta)^2}{n^{-1}\sum_{i=1}^n (Z_i - \theta)^2}
= \chi^2_{1,\;\texttt{conf\_level}}.
$$
Ranking costs
\(O\!\bigl(p\,n\log n\bigr)\); forming and normalising
\(R^\top R\) costs \(O\!\bigl(n p^2\bigr)\) with \(O(p^2)\) additional
memory. The optional jackknife Euclidean-likelihood confidence intervals add
per-pair delete-one recomputation work and are intended for inference rather
than raw-matrix throughput.