“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Can distance-covariance geometry sharpen PLOM in high-noise settings?

Questions

The question. PLOM uses a Gaussian kernel on ambient distances to learn a data manifold. The kernel treats every coordinate the same. Distance covariance detects non-linear dependence between coordinates (the kind that Pearson correlation misses). Can you use dCov among the ambient coordinates to weight the kernel so PLOM ignores the dimensions that are independent of the manifold structure?

What started this

The Wikipedia visual for distance correlation has a grid of scatter plots where Pearson is zero but the two variables are clearly dependent — a circle, a parabola, a sine wave, etc. dCov is exactly the statistic that picks those up. PLOM works by building a kernel out of pairwise distances and reading off the manifold from its eigenstructure. The two methods both care about joint geometric structure between variables; pairing them should mean something. The naive version: weight PLOM's Gaussian kernel so dimensions that are non-trivially tied together in the dCov sense get more say than dimensions that look like noise to the rest of the data.

Setup: PLOM and dCov in one paragraph each

PLOM (Soize & Ghanem, 2016) generates samples on a learned data manifold. Build a Gaussian kernel on the training data, row-normalize, eigendecompose to get a diffusion-map embedding. Top non-trivial eigenvectors are intrinsic manifold coordinates. New samples are produced by random walks (in the original paper, an Itô SDE; in the implementation on this site, a simpler chord-interpolation scheme) that respect the learned geometry.

Distance covariance (Székely, Rizzo, Bakirov, 2007) is a model-free dependence measure. Compute pairwise distance matrices on each variable, double-center each, take the elementwise product, average. The square root of that average is . The key property: iff and are independent — regardless of linearity. The normalized version, distance correlation, lives in .

The pairing. For data in , you can compute a matrix where entry is . That matrix describes the dependence graph between ambient coordinates — which ones are tied together by the manifold, which ones aren't. The hybrid idea: use this dependence graph to weight PLOM's kernel so that coordinates the rest of the data don't see (the "noise dimensions") get downweighted in the distance.

Test setup

A 2D noisy unit circle in , then embedded in by padding with three pure Gaussian-noise dimensions. The first two coordinates carry the manifold; the other three carry no information about it. samples, circle noise , noise-dimension std .

The test: how tightly do new samples produced by PLOM hug the unit circle in the first two coordinates (where "rms distance from circle" = the radial mean-squared error against )? Training data has rms distance 0.048 — that's the actual noise scale of the circle. A good sampler should match this. Vanilla PLOM, which doesn't know which dimensions matter, has to find the manifold in the noisy 5D ambient space.

Three variants

  1. Vanilla. Standardize each dimension to unit variance. Build the Gaussian kernel on raw 5D distances. Baseline.
  2. Diagonal-weighted. Compute the dCor matrix. For each dimension , define an "importance" score as the max dCor with any other dimension. Square it, normalize so the max is 1, and use the resulting per-dim weights to scale the standardized coordinates before building the kernel.
  3. Eigenbasis. Compute the dCor matrix, subtract the identity (this strips out the uninformative diagonal of 1's), eigendecompose. Project the data into the top positive-eigenvalue eigenvectors with sqrt-eigenvalue scaling. This is the cleanest version of "rotate into the dependence-structure basis."

The numbers

"""
Hybrid dCov-weighted PLOM. Test: 2D circle embedded in 5D, three of those
dimensions pure Gaussian noise.

Three variants:
  A. vanilla PLOM             — Gaussian kernel on standardized X
  B. diagonal-weighted PLOM   — scale each X-dim by its max-dCor importance
  C. eigenbasis PLOM          — eigendecomp of (dCor - I), project onto top-k

Key object: distance correlation MINUS identity. The raw dCor matrix has 1's
on the diagonal (uninformative); subtracting identity isolates pairwise
dependence structure.
"""
import numpy as np

# --- data ----------------------------------------------------------------
rng = np.random.default_rng(0)
N        = 300
NOISE    = 0.05
NOISE_AMBIENT = 0.5     # std of the three garbage dimensions
theta = rng.uniform(0.0, 2 * np.pi, N)
x = np.cos(theta) + NOISE * rng.standard_normal(N)
y = np.sin(theta) + NOISE * rng.standard_normal(N)
noise_dims = rng.standard_normal((N, 3)) * NOISE_AMBIENT
X = np.column_stack([x, y, noise_dims])
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

# --- distance covariance / correlation -----------------------------------
def dcov(u, v):
    A = np.abs(u[:, None] - u[None, :])
    B = np.abs(v[:, None] - v[None, :])
    A = A - A.mean(axis=0)[None, :] - A.mean(axis=1)[:, None] + A.mean()
    B = B - B.mean(axis=0)[None, :] - B.mean(axis=1)[:, None] + B.mean()
    val = np.mean(A * B)
    return np.sqrt(val) if val > 0 else 0.0

def dcor_matrix(X):
    d = X.shape[1]
    dvars = np.array([dcov(X[:, j], X[:, j]) for j in range(d)])
    R = np.eye(d)
    for j in range(d):
        for k in range(j + 1, d):
            denom = np.sqrt(dvars[j] * dvars[k]) + 1e-12
            R[j, k] = dcov(X[:, j], X[:, k]) / denom
            R[k, j] = R[j, k]
    return R

R = dcor_matrix(X_std)

# --- variant B: diagonal weighting by max-dCor importance ----------------
R_off = R.copy(); np.fill_diagonal(R_off, 0.0)
importance = R_off.max(axis=1)
weights = importance ** 2; weights /= weights.max()
X_B = X_std * weights[None, :]

# --- variant C: eigenbasis projection of (dCor - I) ----------------------
lam, U = np.linalg.eigh(R - np.eye(R.shape[0]))
idx = np.argsort(lam)[::-1]
lam, U = lam[idx], U[:, idx]
keep = lam > 1e-3
X_C = X_std @ U[:, keep] * np.sqrt(lam[keep])[None, :]

# --- PLOM diffusion-map + chord-interpolation sampler --------------------
def pairwise_sq(A, B):
    return (np.sum(A**2, axis=1)[:, None]
          + np.sum(B**2, axis=1)[None, :] - 2.0 * A @ B.T)

def diffusion_map(X, sigma, m_eigs):
    K = np.exp(-pairwise_sq(X, X) / (2.0 * sigma ** 2))
    d = K.sum(axis=1); Dsq_inv = 1.0 / np.sqrt(d)
    S = Dsq_inv[:, None] * K * Dsq_inv[None, :]
    w, V = np.linalg.eigh(S)
    idx = np.argsort(w)[::-1]
    return (Dsq_inv[:, None] * V[:, idx])[:, 1:1 + m_eigs]

def manifold_sample(X, phi, n_new, k_neighbors, alpha_range, rng):
    new = np.empty((n_new, X.shape[1]))
    for i in range(n_new):
        i0 = rng.integers(len(X))
        d2 = np.sum((phi - phi[i0]) ** 2, axis=1)
        knn = np.argsort(d2)[1:k_neighbors + 1]
        i1 = knn[rng.integers(k_neighbors)]
        a = rng.uniform(-alpha_range, alpha_range)
        new[i] = X[i0] + a * (X[i1] - X[i0])
    return new

def median_pair(Z):
    sq = pairwise_sq(Z, Z)
    return float(np.sqrt(np.median(sq[np.triu_indices(len(Z), 1)])))

phi_A = diffusion_map(X_std, sigma=median_pair(X_std) * 0.5, m_eigs=4)
phi_B = diffusion_map(X_B,   sigma=median_pair(X_B)   * 0.5, m_eigs=4)
phi_C = diffusion_map(X_C,   sigma=median_pair(X_C)   * 0.5, m_eigs=4)

n_new = 500
gen_A = manifold_sample(X, phi_A, n_new=n_new, k_neighbors=6, alpha_range=1.0, rng=rng)
gen_B = manifold_sample(X, phi_B, n_new=n_new, k_neighbors=6, alpha_range=1.0, rng=rng)
gen_C = manifold_sample(X, phi_C, n_new=n_new, k_neighbors=6, alpha_range=1.0, rng=rng)

# --- measure on-manifold-ness in the 2D signal subspace ------------------
def rms_circle(P, R=1.0):
    r = np.sqrt(P[:, 0]**2 + P[:, 1]**2)
    return np.sqrt(np.mean((r - R) ** 2))

for name, P in [("training",   X),
                ("A vanilla",  gen_A),
                ("B diag-wt",  gen_B),
                ("C eigen",    gen_C)]:
    print(f"{name:<12} rms circle dist = {rms_circle(P):.4f}")
Distance-correlation matrix (1 on diagonal, dCor in [0,1] off-diag):
  +1.000  +0.208  +0.080  +0.079  +0.103
  +0.208  +1.000  +0.093  +0.113  +0.075
  +0.080  +0.093  +1.000  +0.088  +0.094
  +0.079  +0.113  +0.088  +1.000  +0.097
  +0.103  +0.075  +0.094  +0.097  +1.000

Per-dim max-dCor importance: [0.208 0.208 0.094 0.113 0.103]
Per-dim weights (importance^2, normalized): [1.000 1.000 0.203 0.298 0.244]

Eigenvalues of (dCor - I), descending: [+0.418 -0.024 -0.087 -0.091 -0.217]
Keeping 1 positive-eigenvalue components.

variant                                rms dist    max dist   noise std
------------------------------------------------------------------------
training data                            0.0480      0.1619      0.5046
A. vanilla PLOM (standardized)           0.0878      0.6884      0.5795
B. diagonal-weighted PLOM                0.0579      0.2396      0.6621
C. eigenbasis PLOM (dCor - I)            0.4672      1.9635      0.6336

Three things to read off. The dCor matrix correctly picks out the signal: entry is 0.208, all other off-diagonal entries sit in the range 0.075–0.113. The signal-to-noise ratio in dCor space is about 2×.

Variant B wins. RMS distance from the circle: vanilla 0.088, diagonal-weighted 0.058 — a 35% reduction. Diagonal weighting is the simplest possible use of the dCor matrix and it works. The per-dim importance scores are 0.21 for the two signal coordinates, 0.09–0.11 for the three noise coordinates. After squaring and normalizing, the signal coordinates keep weight 1.0 and the noise coordinates drop to 0.20–0.30. The kernel still sees all five dimensions, but distances along the noise dimensions contribute less.

Variant C — the one closest to the "use the dCov matrix as a kernel" framing I started with — failed. RMS 0.467. The reason: has only one positive eigenvalue at this sample size (0.42). The other four are slightly negative — sampling artifacts, not real "anti-dependence" structure. Projecting the data into a 1D subspace before running PLOM threw away too much for the chord-interpolation sampler to recover anything sensible. The eigenbasis approach only makes sense when there are multiple positive eigenvalues — i.e., when there are multiple genuinely-dependent groups of coordinates, not just one signal pair against a noise floor.

What this tells me

Diagonal weighting works because it's a soft mask — every coordinate still contributes to the kernel, just by different amounts. The signal coordinates get the full weight they'd have in vanilla PLOM; the noise coordinates get a quarter of it. The chord-interpolation sampler in ambient space still has full information, just with a manifold-detection geometry that favors the right subspace. Eigenbasis projection is the wrong move when only one dependence direction is detectable, because it hard-truncates information the sampler still needs.

The dCor signal also has to be detectable at the sample sizes you're using. At with these parameters the signal-to-noise contrast in dCor space is about 2×. I tried an earlier version at and the contrast dropped to ~1.4×; neither variant did better than vanilla in that regime. So the hybrid lives in the middle: signal detectable above the finite-N dCor floor, but not so dominant that vanilla PLOM trivially recovers it.

The structural finding worth keeping is that is the right object for any future use of dCov as a kernel between coordinates. The raw dCov matrix is dominated by self-dCov on the diagonal, which only measures per-dimension spread. The dCor matrix has 1's on the diagonal which are uninformative. Subtracting identity isolates exactly the off-diagonal pairwise dependence structure.

Open threads

  1. Swiss roll. A real test for the eigenbasis variant. The Swiss roll has multi-scale dependence — distant arms of the spiral are correlated via the underlying parameter, and Euclidean distance across the gap is misleading. dCov should pick this up better than Gaussian distance, and there might be more than one positive eigenvalue of .
  2. Sample-size scaling. The hybrid's win at is 35%. Does the gap grow with (better dCor estimate → sharper weights) or shrink (vanilla PLOM has more data to work with anyway)?
  3. The Brownian-distance kernel connection. Sejdinovic, Sriperumbudur, Gretton, Fukumizu (2013) showed dCov is equivalent to HSIC with a specific characteristic kernel (the "Brownian-distance kernel"). The diagonal-weighted scheme here is approximating a particular use of that kernel — there might be a cleaner formulation that drops out of writing the whole thing in the HSIC framework rather than building diagnostics by hand.
  4. Real high-dimensional data. The test here is synthetic. Whether the hybrid helps on real data — turbulence fields, molecular conformations, audio — depends on whether real datasets have the "many noise dimensions, few signal dimensions" structure the test was built for.

References I want to chase down