Local Correlation Coefficients on Graphs

Compute vertex-level correlation coefficients measuring the alignment of directional changes between functions defined on graph vertices. This function provides a unified interface supporting vector-vector, vector-matrix, and matrix-matrix computations with automatic dispatch based on input types.

Usage

lcor(
  adj.list,
  weight.list,
  y,
  z,
  type = c("derivative", "unit", "sign"),
  y.diff.type = c("difference", "logratio"),
  z.diff.type = c("difference", "logratio"),
  epsilon = 0,
  winsorize.quantile = 0,
  instrumented = FALSE,
  mc.cores = 1L,
  hop.radius = 1L,
  pairs = NULL
)

Arguments

adj.list

List of integer vectors containing 1-based vertex indices. Element i contains the neighbors of vertex i.

weight.list

List of numeric vectors containing edge weights (lengths). Must have same structure as adj.list.

y

Numeric vector, matrix, or data frame of function values. For vectors, length must equal the number of vertices. For matrices/data frames, number of rows must equal the number of vertices.

z

Numeric vector, matrix, or data frame of function values. Dimension requirements match those for y.

type

Character scalar specifying the weighting scheme:

"derivative": Geometric weights (w_e = 1/length^2). Normalizes by edge length squared, making the measure analogous to comparing derivatives. Appropriate when functions represent continuous quantities sampled at irregular spatial positions. This is the default.
"unit": Equal weights (w_e = 1). Treats all edges equally regardless of length, while retaining the magnitudes of the field differences.
"sign": Legacy name currently equivalent to "unit". It retains edge-difference magnitudes; it is not a sign-only statistic.

y.diff.type

Character scalar specifying edge difference type for y:

"difference": Standard differences: Delta_e f = f(u) - f(v). Appropriate for continuous data in Euclidean space.
"logratio": Log-ratios: Delta_e f = log((f(u) + epsilon) / (f(v) + epsilon)). A per-edge transform for nonnegative fields; see its limitations below. The log transformation maps multiplicative changes to an additive scale.

z.diff.type

Character scalar specifying edge difference type for z. Same options as y.diff.type.

epsilon

Numeric scalar for pseudocount in log-ratio transformations. Must be finite and nonnegative. Zero selects 1e-6 times the smallest positive value (1e-6 when all values are zero). Vector-vector and matrix-matrix paths select it separately per field; vector-matrix selects one for y and one across the whole feature matrix. Supply an explicit positive value when comparing dispatch paths or chunks using log-ratios. Only used for log-ratios.

winsorize.quantile

Numeric scalar for winsorization of edge differences. If 0 (default), no winsorization is applied. If positive (e.g., 0.025), clips edge differences to the $[q, 1-q]$ percentile range for robustness against outliers.

instrumented

Logical scalar. If TRUE and both y and z are vectors, returns additional diagnostic information including edge-level differences and weights. Ignored for matrix inputs. Default is FALSE.

mc.cores

Integer scalar specifying the number of cores for parallel computation in the matrix-matrix case. Default is 1 (sequential). Note that parallel::mclapply is not available on Windows.

hop.radius

Integer scalar specifying the hop distance for local neighborhoods. Default is 1 (immediate neighbors). When hop.radius > 1, neighborhoods include all vertices reachable within hop.radius hops. For type = "derivative", edge lengths are replaced by shortest path lengths (sum of edge lengths) among paths with at most hop.radius hops. For type = "unit" or type = "sign", all k-hop edges are treated as having length 1.

pairs

Optional two-column matrix of column indices for matrix y and z. Computes only these pairs, in supplied order, as a vertex-by-pair matrix. Repeated pairs and self-pairs are allowed. Column names identify both fields; the pairs attribute records the exact indices. Uses serial vector-pair computation, regardless of mc.cores; requires instrumented = FALSE.

Value

The return type depends on the input types:

Vector-vector (y and z both vectors): An object of class "lcor_vector_vector_result". If instrumented = FALSE, a numeric vector of correlation coefficients at each vertex, with values in $[-1, 1]$. If instrumented = TRUE, a list containing:

vertex.coefficients: Numeric vector of local correlation coefficients.
vertex.delta.y: List of edge differences for y at each vertex.
vertex.delta.z: List of edge differences for z at each vertex.
vertex.weights: List of edge weights at each vertex.
all.delta.y: All y edge differences (pre-winsorization).
all.delta.z: All z edge differences (pre-winsorization).
y.lower, y.upper: Winsorization bounds for y.
z.lower, z.upper: Winsorization bounds for z.

Vector-matrix (y vector, z matrix) or matrix-vector (y matrix, z vector): An object of class "lcor_vector_matrix_result". If instrumented = FALSE, a numeric matrix of dimension (n.vertices x n.columns) where column j contains local correlation coefficients between the vector and the j-th column of the matrix. If instrumented = TRUE, a list containing:

column.coefficients: Matrix (n.vertices x n.columns) of local correlations.
y.lower, y.upper: Scalar winsorization bounds for the vector input.
z.lower, z.upper: Numeric vectors of length n.columns giving per-column winsorization bounds for the matrix input.

For matrix-vector input, the attribute "transposed" is set to TRUE.

Matrix-matrix with identical(y, z) = TRUE (symmetric case): An object of class "lcor_matrix_matrix_result" with attribute symmetric = TRUE. Only upper triangular pairs (i, j) with i < j are computed. If instrumented = FALSE, a numeric matrix of dimension (n.vertices x n.pairs) where n.pairs = ncol(y) * (ncol(y) - 1) / 2. Column names are in "col_i:col_j" format. If instrumented = TRUE, a list containing:

pair.coefficients: Matrix (n.vertices x n.pairs) of local correlations.
pair.names: Character vector of pair names in "col_i:col_j" format.
pair.indices: Integer matrix (n.pairs x 2) of column index pairs.
column.lower, column.upper: Numeric vectors of winsorization bounds for each column (shared between y and z since they are identical).

Matrix-matrix with identical(y, z) = FALSE (asymmetric case): An object of class "lcor_matrix_matrix_result" with attribute symmetric = FALSE. All pairs (i, j) for i in 1:ncol(y) and j in 1:ncol(z) are computed. If instrumented = FALSE, a 3D numeric array of dimension (n.vertices x ncol(y) x ncol(z)). Element [v, i, j] contains lcor(y[,i], z[,j]) at vertex v. If instrumented = TRUE, a list containing:

pair.coefficients: 3D array (n.vertices x ncol(y) x ncol(z)) of local correlations.
y.column.lower, y.column.upper: Numeric vectors of winsorization bounds for each column of y.
z.column.lower, z.column.upper: Numeric vectors of winsorization bounds for each column of z.

Details

The local correlation coefficient at vertex v measures the alignment of directional changes in y and z within v's neighborhood:

$$lcor(y,z)(v) = \frac{\sum_e w_e \Delta_e y \cdot \Delta_e z} {\sqrt{\sum_e w_e (\Delta_e y)^2} \sqrt{\sum_e w_e (\Delta_e z)^2}}$$

where the sum is over edges e incident to vertex v, and Delta_e f denotes the edge difference (either standard or log-ratio). This is an uncentered weighted cosine of incident edge differences, not Pearson correlation across vertices. Isolates and fields with weighted local norm at most 1e-10 return zero; that convention is not evidence of orthogonality. Derivative weighting omits edges with lengths at most 1e-10. A nonzero epsilon and clipping can alter rescaling properties. Local alignment alone does not establish causation.

Dispatch Behavior

The function automatically dispatches to specialized implementations based on input types:

Case (a) Vector-Vector: Computes a single local correlation between two functions. This is the fundamental operation underlying all other cases.

Case (b) Vector-Matrix: Efficiently computes correlations between one vector y and each column of matrix Z. The implementation pre-computes y-dependent quantities once and reuses them across columns, providing significant speedup over repeated vector-vector calls.

Case (c) Matrix-Vector: Equivalent to case (b) with roles swapped. The function internally calls the vector-matrix implementation with y and z exchanged and diff.types appropriately swapped.

Case (d) Matrix-Matrix: Computes correlations between all pairs of columns from y and z. When y and z are identical (checked via identical()), only upper triangular pairs are computed to avoid redundancy. Parallel computation via mc.cores can accelerate this case substantially.

Edge Difference Types

Ordinary differences describe absolute changes in a field across an edge. Log-ratios describe multiplicative changes after adding the pseudocount. Either field can use either transform if its input domain permits it. Choosing a transform is a modeling assumption; relative-abundance data alone do not establish that edge log-ratios supply an appropriate compositional analysis.

Geometric Interpretation

With ordinary differences and no winsorization, for linear fields with gradients a and b on an embedded star, define $M = \sum_e w_e d_e d_e^T$, where d_e is its displacement vector. Substituting $\Delta_e y = a^T d_e$ in the discrete formula gives $a^T M b / \sqrt{(a^T M a)(b^T M b)}$. This equals the Euclidean gradient cosine when M is proportional to the identity on the relevant span. A smooth-field limit additionally needs shrinking neighborhoods, controlled Taylor remainders, nonzero gradients, and that isotropic directional sampling. Shrinking an anisotropic star alone does not establish this limit.

For example, neighbors (1,0) and (1,1) around (0,0), with fields x and y, give local alignment 1/sqrt(3) at the center under derivative weighting, although the Euclidean gradients are orthogonal. This follows directly from the displayed finite-graph formula.

A log-ratio here compares one field between adjacent vertices after adding the stated pseudocount; it does not compare different components within a composition. Zeros and the pseudocount affect the result, and this transform alone is not a general compositional-inference guarantee.

Hop Radius

By default, the local correlation at vertex v is computed using edges incident to v (hop.radius = 1). For hop.radius > 1, the neighborhood expands to all vertices within hop.radius hops of v. The computation is then performed on the k-hop star graph centered at v. For type = "derivative", k-hop edges are weighted by the inverse squared shortest path length within the hop limit; for type = "unit" and type = "sign", all k-hop edges are treated as having unit length.

Memory and bounded calculations

With n vertices and p and q fields, numeric output alone needs approximately 8n bytes for vectors, 8np for vector-matrix, 8np(p-1)/2 for identical matrix inputs, and 8npq for distinct matrices. Inputs, expanded hop graphs, worker copies, pair-result lists, and assembly temporaries add to these totals. For example, n=10000 and p=q=100 needs 800 MB for the full asymmetric output. pairs reduces output to 8n*K bytes for K requested pairs, plus one temporary vertex vector. To bound a long analysis, split the requested pairs into small blocks, call lcor(..., pairs = block), save each returned block, then discard it before the next call. Save the pair indices with every block and record completed blocks; a stopped analysis must not be treated as complete. Use an explicit common positive epsilon when comparing logratio results across dispatch paths or chunks. Native kernels do not currently poll for interrupts inside a single pair; R can process interrupts between pair calls.

References

Pawel Gajer and Jacques Ravel (2025). From Global to Local Correlation: Geometric Decomposition of Statistical Inference. arXiv:2511.04599.

Examples

# Three vertices, with one field increasing in the same direction as y.
adjacency <- list(c(2L, 3L), c(1L, 3L), c(1L, 2L))
lengths <- lapply(adjacency, function(v) rep(1, length(v)))
y <- c(0, 1, 3)
features <- cbind(aligned = 2 * y + 1, opposed = -y)

# Default vector output: +1 means aligned incident differences.
alignment <- lcor(adjacency, lengths, y, features[, "aligned"])
as.numeric(alignment)
#> [1] 1 1 1
stopifnot(isTRUE(all.equal(as.numeric(alignment), rep(1, 3))))

# Matrix output has vertices in rows and features in columns.
coefficients <- lcor(adjacency, lengths, y, features)
colMeans(coefficients) # aligned = +1; opposed = -1
#> aligned opposed 
#>       1      -1 
stopifnot(identical(dim(coefficients), c(3L, 2L)))

# Request the diagnostic list explicitly when inspecting edge quantities.
diagnostic <- lcor(adjacency, lengths, y, features[, "aligned"],
                   instrumented = TRUE)
diagnostic$vertex.coefficients
#> [1] 1 1 1
diagnostic$vertex.weights
#> [[1]]
#> [1] 1 1
#> 
#> [[2]]
#> [1] 1 1
#> 
#> [[3]]
#> [1] 1 1
#> 

# Identical matrices return one column for each distinct feature pair.
pairs <- lcor(adjacency, lengths, features, features)
colnames(pairs)
#> [1] "aligned:opposed"
# These are descriptive local alignments, not significance tests.