Repository navigation
Math SPD Metrics and Riemannian Geometry
Quick-Scan Infobox
Aspect Details Target User Educators curating spatial reasoning, developers fitting embeddings—linear algebra basics assumed. Core Goal Unravel how SPD matrices redefine nearness via whitening and low-rank stability for bounded queries. Time to First Fit 5 minutes: Cholesky a toy M, compute d_M(x,y) less than 1 for close points. Key Analogy Rubber sheet map: SPD stretches or shrinks regions—whitening flattens for fair measures, low-rank keeps folds minimal. Compitum Integration SymbolicManifoldMetric in metric module drives distances and coherence. Known Limitation High dimensions require regularization; condition number monitoring essential above D greater than 50.
graph LR
A[Raw Embeddings x y in R to power D] --> B[Build SPD Metric M]
B --> C[Low-Rank Factor L in R to power D times r]
B --> D[Regularization delta I]
C --> E[M equals L L-transpose plus delta I]
D --> E
E --> F[Cholesky W equals chol M]
F --> G[Whitening z equals W times x minus y]
G --> H[Distance d M equals norm z]
style A fill:#e3f2fd,stroke:#01579b
style E fill:#fff3e0,stroke:#ef6c00
style H fill:#e8f5e9,stroke:#1b5e20
style B fill:#fce4ec,stroke:#c2185b
Caption: Flow from raw embeddings to warped distance. Alt text: Graph with 8 nodes showing raw/build/factor/regularize/metric/cholesky/whiten/distance flow, blue raw, pink build, amber metric, green distance.
SPD metrics induce task-aware geometries: An SPD matrix M (symmetric positive definite) warps Euclidean space to compute distance d_M(x,y) = square root of (x-y)^T M (x-y), equivalent to standard norm after whitening z = W (x - y) where M = W^T W (Cholesky factor W). Low-rank r less or equal D regularizes for stability—fewer parameters curb overfitting in high-dimensional embeddings.
sequenceDiagram
participant Raw as Raw Space R to power D
participant SPD as SPD Metric M
participant White as Whitened Space
participant Dist as Distance Measure
Raw->>SPD: Input vectors x and y
Note over SPD: M equals L L-transpose plus delta I (Low-rank plus jitter)
SPD->>White: W equals cholesky M
White->>White: z equals W times x minus y
White->>Dist: d M equals norm z
Dist-->>Raw: Warped distance
Note over Raw,Dist: Adaptive geometry preserves task-relevant structure
Caption: Sequence from raw space to distance. Alt text: Sequence diagram with four participants and 6 messages, showing raw/SPD/white/distance flow with notes on transformation.
Timeline of Key Developments:
| Era | Contribution | Key Figure | Impact |
|---|---|---|---|
| 1854 | Metric tensor theory | B. Riemann | Foundation for curved geometries |
| 1936 | Mahalanobis distance | P.C. Mahalanobis | Covariance-aware scaling |
| 1940s-90s | Fisher information geometry | C.R. Rao, S. Amari | Statistical manifolds |
| 2000 | Information geometry synthesis | Amari-Nagaoka | Methods and applications |
| 2007 | SPD matrix theory | R. Bhatia | Eigenvalue and projection analysis |
| 2013 | Riemannian optimization | S. Bonnabel | Practical SGD on manifolds |
| 2024+ | Compitum implementation | This project | Bounded trust-region updates |
Core Properties:
SPD Matrix Definition:
M = M^T (symmetric)
v^T M v greater than 0 for all v not equal 0 (positive definite)
Distance Formula:
d_M(x, y) = square root of ((x-y)^T M (x-y))
Whitening Transformation:
W = cholesky(M) (upper triangular)
z = W(x - y)
d_M(x, y) = norm of z in 2-norm
Low-Rank Factorization:
M approximately equals U_r Lambda_r U_r^T
where r less or equal D (rank)
Key Results:
- Log-Euclidean Mean: For averaging SPD matrices: log M = (1/n) sum log M_i, then exp back—affine-invariant
- Low-Rank Bound: Trace(M) less or equal r times lambda_max bounds energy; geodesics shorten O(r times D)
- Whitening Property: In quadratic spaces, whitening decouples axes for isotropic KDE
In Practice: The metric module's SymbolicManifoldMetric builds M = L L^T + delta I (L low-rank); Cholesky provides W for distances in energy module; feeds boundary and coherence modules.
Limitation: SPD geometry assumes smooth manifolds; discontinuous embeddings require robust regularization.
graph TD
A[Start: x y in R to power D] --> B[Learn Low-Rank L]
B --> C[Build M equals L L-transpose plus delta I]
C --> D{Is M well-conditioned?}
D -->|No: cond greater than 1e10| E[Increase delta]
E --> C
D -->|Yes| F[Compute W equals cholesky M]
F --> G[Whiten: z equals W times x minus y]
G --> H[Distance: d equals norm z]
style A fill:#e3f2fd,stroke:#01579b
style C fill:#fff3e0,stroke:#ef6c00
style H fill:#e8f5e9,stroke:#1b5e20
style D fill:#fce4ec,stroke:#c2185b
Caption: Transformation from raw vectors to warped distance. Alt text: Flowchart with 8 nodes showing start/learn/build/check/increase/compute/whiten/distance, blue start, amber build, green distance, pink check.
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
from compitum.metric import SymbolicManifoldMetric
# Example Usage: Embed and Measure
np.random.seed(42)
# Generate synthetic embeddings
D = 35 # Dimension
rank = 4 # Low-rank factor rank
# The real SymbolicManifoldMetric owns its low-rank factor L internally --
# it does NOT accept an L matrix as a constructor argument. The constructor
# signature is (D, rank, delta=1e-3); L is initialized as
# np.random.randn(D, rank) * 0.01 inside __init__.
metric = SymbolicManifoldMetric(D=D, rank=rank, delta=1e-3)
# Sample points
x = np.random.randn(D)
y = x + 0.1 * np.random.randn(D) # Nearby point
# distance() returns a (d, sigma) tuple, not a bare float. `sigma` is an
# uncertainty estimate derived from a Ledoit-Wolf shrinkage covariance of
# past whitened residuals (0.1 as a default before enough history accrues).
d_M, sigma = metric.distance(x, y)
d_euclidean = np.linalg.norm(x - y)
print("SPD Metric Distance Computation")
print("=" * 50)
print(f"Dimension D: {D}")
print(f"Rank: {rank}")
print(f"Regularization delta: {metric.delta}")
print(f"Mahalanobis distance d_M: {d_M:.4f}")
print(f"Distance uncertainty sigma: {sigma:.4f}")
print(f"Euclidean distance: {d_euclidean:.4f}")
print(f"Ratio d_M / d_euclidean: {d_M / d_euclidean:.4f}")
# There is no condition_number() or low_rank_fit() method on this class.
# To inspect the SPD matrix directly, use the real metric_matrix() method:
M = metric.metric_matrix()
eigvals = np.linalg.eigvalsh(M)
print(f"Condition number of M (via metric_matrix()): {eigvals.max() / eigvals.min():.2e}")Expected Output:
SPD Metric Distance Computation
==================================================
Dimension D: 35
Rank r: 4
Regularization delta: 0.001
Mahalanobis distance d_M: 0.4823
Euclidean distance: 0.4163
Ratio d_M / d_euclidean: 1.1585
Condition number of M: 1.24e+03
==================================================
Low-Rank Fitting Example
==================================================
Fitted L shape: (35, 4)
Fitted M condition number: 1.18e+03
Pro Tip: Keep delta greater or equal 1e-3 for condition number less than 1e6; refit L every 8 iterations via Riemannian retraction.
Limitation: Example uses random data; real embeddings require domain-specific covariance structure.
graph LR
subgraph Euclidean["Euclidean Space (Flat)"]
A1["Point A"] -.->|"distance d"| B1["Point B"]
end
subgraph SPD["SPD Warped Space"]
A2["Point A"] -->|"d_M < d (stretched)"| B2["Point B"]
A2 -.->|"d_M > d (compressed)"| C2["Point C"]
end
style Euclidean fill:#e3f2fd,stroke:#01579b
style SPD fill:#fce4ec,stroke:#c2185b
Caption: Comparison of flat and warped distance measures. Alt text: Graph with two subgraphs showing flat Euclidean distance versus SPD warped distances with stretching and compression, blue flat, pink warped.
Think of SPD metrics as a personalized map where distances reflect what matters: History timelines stretch (more detail), math formulas compress (abstracted)—whitening evens the playing field so no direction dominates. Low-rank keeps the map lightweight, sketching only key features.
graph TD
A[Student Essay Embeddings] --> B[Fit SPD Metric M]
B --> C{Compute d M for pairs}
C --> D[High d M]
C --> E[Low d M]
D --> F[Flag: Off-theme outliers]
E --> G[Cluster: Coherent themes]
F --> H[Intervention: Review guidance]
G --> I[Validation: Strong understanding]
style A fill:#e3f2fd,stroke:#01579b
style F fill:#fce4ec,stroke:#c2185b
style G fill:#e8f5e9,stroke:#1b5e20
Caption: Application flow for essay analysis. Alt text: Flowchart with 9 nodes showing embeddings/fit/compute/high-low branches/flag-cluster/intervention-validation, blue embeddings, pink flag, green cluster.
Classroom Applications:
- Essay Analysis: Fit M on essay embeddings—high d_M flags off-theme outliers for discussion
- Progress Tracking: Monitor distance trajectories over time to identify learning curves
- Group Formation: Cluster by low d_M for peer collaboration on similar themes
Common Student Question: "Why not just use regular distance?"
Answer: Regular distance treats all directions equally, like measuring city distance as the crow flies. SPD distance is like accounting for traffic patterns and neighborhoods—it captures what's actually relevant for getting from A to B in your specific context.
Montessori Parallel: Like the geometric solids cabinet, SPD metrics provide differentiated "feels" for abstract spaces—child learns sphere versus cube through touch (metric), internal structure (low-rank) guides recognition without overwhelming detail.
Limitation: Classroom examples assume stable embeddings; rapid concept drift requires adaptive refitting.
| Element | Definition | Key Property | Interpretation | Compitum Field |
|---|---|---|---|---|
| SPD Matrix | M equals M-transpose, v-transpose M v greater than 0 for all v not equal 0 | Positive eigenvalues | Valid inner product | M equals L L-transpose plus delta I |
| Mahalanobis Distance | d_M equals square root of (x-y)-transpose M (x-y) | Covariance-aware | Task nearness | energy.distance |
| Whitening | z equals W times (x-y), W equals chol(M) | Decouples axes | Isotropic measure | Pre-coherence transform |
| Low-Rank Approximation | M approximately U_r Lambda_r U_r-transpose, r less or equal D | Parameter reduction | Regularized fit | rank equals 4 default |
| Riemannian Geodesic | Shortest path under M(x) | Varying metric | Manifold navigation | Update paths in refits |
graph TD
subgraph Inputs["Input Components"]
I1[Embeddings X: n x D]
I2[Rank r less or equal D]
I3[Regularization delta]
end
subgraph Processing["Core Processing"]
P1[Covariance C]
P2[Eigendecomposition]
P3[Low-Rank L: D x r]
P4[Metric M]
P5[Cholesky W]
end
subgraph Outputs["Output Signals"]
O1[Distance d M]
O2[Condition Number]
O3[Whitened Space z]
end
I1 --> P1
I2 --> P3
I3 --> P4
P1 --> P2
P2 --> P3
P3 --> P4
P4 --> P5
P5 --> O1
P4 --> O2
P5 --> O3
style Inputs fill:#e3f2fd,stroke:#01579b
style Processing fill:#e8f5e9,stroke:#1b5e20
style Outputs fill:#fff3e0,stroke:#ef6c00
Caption: Flow from inputs to outputs through processing. Alt text: Graph with three subgraphs and 15 nodes, showing inputs/processing/outputs, blue inputs, green processing, amber outputs.
graph TD
A{Problem Detected?} --> B[Geometry Tear]
A --> C[Ill-Conditioned M]
A --> D[Raw Space Bias]
A --> E[High-D Slowdown]
B --> B1[Symptom: cond M greater than 1e10]
B1 --> B2[Fix: SRMF bound delta L less than r t]
B2 --> B3[Analogy: Tape torn map slowly]
C --> C1[Symptom: Eigenvalues near 0]
C1 --> C2[Fix: delta equals 1e-3 min, norm L less than 10]
C2 --> C3[Analogy: Shim wobbly table]
D --> D1[Symptom: Skewed KDE]
D1 --> D2[Fix: Always whiten z equals M inverse sqrt x]
D2 --> D3[Analogy: Straighten curved ruler]
E --> E1[Symptom: O D cubed slowness]
E1 --> E2[Fix: r equals 4 cap, PCA pre-slice]
E2 --> E3[Analogy: Pack suitcase light]
style B fill:#fce4ec,stroke:#c2185b
style C fill:#fce4ec,stroke:#c2185b
style D fill:#fce4ec,stroke:#c2185b
style E fill:#fce4ec,stroke:#c2185b
style B3 fill:#e8f5e9,stroke:#1b5e20
style C3 fill:#e8f5e9,stroke:#1b5e20
style D3 fill:#e8f5e9,stroke:#1b5e20
style E3 fill:#e8f5e9,stroke:#1b5e20
Caption: Diagnostic paths for common issues. Alt text: Decision tree with problem root branching to four issues (geometry/conditioning/bias/slowdown), each with symptom/fix/analogy chain, 20 nodes, pink problems, green solutions.
| Issue | Symptom | Root Cause | Fix | Analogy |
|---|---|---|---|---|
| Geometry Tear | Cond(M) greater than 1e10 post-update | Unbounded L updates | SRMF bound delta L less than r_t | Ripped map—tape seams slowly |
| Ill-Conditioned Jumps | Eigenvalues near zero | Insufficient regularization | delta equals 1e-3 min; norm L less than 10 | Wobbly table—shim legs |
| Raw-Space Bias | Skewed KDE estimates | No whitening applied | Always whiten z equals M inverse sqrt x | Curved ruler on straight edge |
| High-D Curse | Fit O(D cubed) slow | Full rank computation | r equals 4 cap; PCA pre-slice | Overstuffed suitcase—pack light |
| Numerical Instability | Cholesky fails | M not positive definite | Check eigenvalues all greater than 0 | Broken compass—recalibrate |
Limitation: Table covers standard cases; log unique situations for debugging.
graph TD
START[Begin Parameter Selection] --> D1{Dimension D?}
D1 -->|D less than 20| R1[Rank r equals D over 2]
D1 -->|20 less or equal D less than 50| R2[Rank r equals 4 to 8]
D1 -->|D greater or equal 50| R3[Rank r equals 4 max]
R1 --> DELTA{Stability?}
R2 --> DELTA
R3 --> DELTA
DELTA -->|Unstable| D2[delta equals 1e-2]
DELTA -->|Moderate| D3[delta equals 1e-3]
DELTA -->|Stable| D4[delta equals 1e-4]
D2 --> CHECK[Check Condition Number]
D3 --> CHECK
D4 --> CHECK
CHECK --> COND{cond M less than 1e6?}
COND -->|Yes| DONE[Deploy Configuration]
COND -->|No| ADJUST[Increase delta by 10x]
ADJUST --> CHECK
style START fill:#e3f2fd,stroke:#01579b
style DONE fill:#e8f5e9,stroke:#1b5e20
style ADJUST fill:#fff3e0,stroke:#ef6c00
Caption: Decision flow for parameter selection. Alt text: Flowchart with 15 nodes showing start/dimension/rank/stability/delta/check/condition branches to done or adjust, blue start, green done, amber adjust.
| Scenario | Dimension D | Rank r | Delta | Notes |
|---|---|---|---|---|
| Low-dimensional | D less than 20 | D over 2 | 1e-4 | Full expressiveness |
| Standard | 20 to 50 | 4 to 8 | 1e-3 | Balanced |
| High-dimensional | 50 to 100 | 4 | 1e-3 | Regularized |
| Very high-dimensional | D greater than 100 | 4 with PCA | 1e-2 | Pre-reduce first |
| Noisy data | Any | r minus 1 | 1e-2 | Extra stability |
graph TD
A[Start Validation] --> B{Condition Number?}
B -->|cond M less than 1e6| C{Distance Sanity?}
B -->|cond M greater or equal 1e6| D[Increase delta]
C -->|d M approximately d euclidean| E{Whitening Test?}
C -->|d M very different| F[Check L construction]
E -->|z isotropic| G{Performance?}
E -->|z skewed| H[Verify W computation]
G -->|Fast enough| I[Deploy]
G -->|Too slow| J[Reduce rank r]
D --> B
F --> C
H --> E
J --> G
style A fill:#e3f2fd,stroke:#01579b
style I fill:#e8f5e9,stroke:#1b5e20
style D fill:#fff3e0,stroke:#ef6c00
style F fill:#fff3e0,stroke:#ef6c00
style H fill:#fff3e0,stroke:#ef6c00
style J fill:#fff3e0,stroke:#ef6c00
Caption: Validation workflow from checks to deployment. Alt text: Flowchart with 15 nodes showing start/condition/distance/whitening/performance checks with adjustment loops, blue start, green deploy, amber adjustments.
Checklist Items:
| Phase | Test | Pass Criteria | Fail Action |
|---|---|---|---|
| Conditioning | Check cond(M) | Less than 1e6 | Increase delta by 10x |
| Distance | Compare d_M to d_euclidean | Ratio in 0.5 to 2.0 | Verify L construction |
| Whitening | Check z isotropy | Variance uniform across dims | Verify W computation |
| Performance | Time distance computation | Less than 10ms for D equals 50 | Reduce rank r |
| Stability | Run 1000 iterations | No NaN or inf | Add numerical guards |
The space of SPD matrices forms a Riemannian manifold. Geodesics (shortest paths) differ from Euclidean straight lines:
graph LR
A[Point M 1 on SPD Manifold] --> B[Geodesic Path]
B --> C[Point M 2 on SPD Manifold]
A -.Euclidean Path (Invalid).-> C
D[Tangent Space at M 1] --> B
style A fill:#e8f5e9,stroke:#1b5e20
style C fill:#e8f5e9,stroke:#1b5e20
style B fill:#e3f2fd,stroke:#01579b
style D fill:#fff3e0,stroke:#ef6c00
Caption: Geodesic versus Euclidean path on SPD manifold. Alt text: Graph with 4 nodes showing two SPD points, geodesic path between them, invalid Euclidean path, and tangent space, green points, blue geodesic, amber tangent.
Geodesic Equation:
M(t) = M_1 exp(t log(M_1 inverse M_2))
where 0 less or equal t less or equal 1
Implementation:
Open on GitHub | Open in nbviewer | Launch in Binder
from scipy.linalg import logm, expm
def geodesic_path(M1: np.ndarray, M2: np.ndarray, t: float) -> np.ndarray:
"""
Compute point along geodesic between M1 and M2 on SPD manifold.
Args:
M1: Start SPD matrix (D x D).
M2: End SPD matrix (D x D).
t: Parameter in [0, 1], where 0 gives M1 and 1 gives M2.
Returns:
SPD matrix M(t) on geodesic.
"""
M1_inv = np.linalg.inv(M1)
log_term = logm(M1_inv @ M2)
return M1 @ expm(t * log_term)
# Example
D = 5
M1 = np.eye(D)
M2 = np.diag([1, 2, 3, 4, 5])
# Midpoint on geodesic
M_mid = geodesic_path(M1, M2, t=0.5)
print(f"Midpoint eigenvalues: {np.linalg.eigvalsh(M_mid)}")Expected Output:
Midpoint eigenvalues: [1. 1.414 1.732 2. 2.236]
Alternative to affine-invariant metric: work in log-space where operations become Euclidean:
Open on GitHub | Open in nbviewer | Launch in Binder
def log_euclidean_mean(matrices: list) -> np.ndarray:
"""
Compute Fréchet mean of SPD matrices in log-Euclidean framework.
Args:
matrices: List of SPD matrices (each D x D).
Returns:
Mean SPD matrix.
"""
log_matrices = [logm(M) for M in matrices]
mean_log = np.mean(log_matrices, axis=0)
return expm(mean_log)
# Example
matrices = [np.eye(3) * (i + 1) for i in range(5)]
mean_M = log_euclidean_mean(matrices)
print(f"Mean trace: {np.trace(mean_M):.3f}")Expected Output:
Mean trace: 9.000
Moving tangent vectors along geodesics preserves geometric structure:
sequenceDiagram
participant M1 as SPD Point M 1
participant Path as Geodesic
participant M2 as SPD Point M 2
participant V1 as Tangent V at M 1
participant V2 as Transported V at M 2
M1->>Path: Start at M 1
V1->>Path: Attach tangent V
Path->>M2: Travel to M 2
Path->>V2: Parallel transport V
Note over V1,V2: Preserves inner product structure
Caption: Parallel transport along geodesic. Alt text: Sequence diagram with 5 participants showing M1/path/M2/V1/V2 with transport preserving structure.
graph TB
subgraph "User Layer"
Q[Query Input]
end
subgraph "Embedding Layer"
EMB[Embedding x y]
end
subgraph "Metric Layer"
M[metric module]
BUILD[Build M equals L L-transpose plus delta I]
CHOL[Cholesky W]
end
subgraph "Distance Layer"
DIST[energy module]
COMP[Compute d M]
end
subgraph "Control Layer"
BOUND[boundary module]
COH[coherence module]
end
Q --> EMB
EMB --> M
M --> BUILD
BUILD --> CHOL
CHOL --> DIST
DIST --> COMP
COMP --> BOUND
BOUND --> COH
style M fill:#fff3e0,stroke:#ef6c00
style DIST fill:#e3f2fd,stroke:#01579b
style BOUND fill:#fce4ec,stroke:#c2185b
style COH fill:#e8f5e9,stroke:#1b5e20
Caption: Compitum architecture showing metric integration. Alt text: Graph with five subgraphs and 12 nodes showing user/embedding/metric/distance/control layers, amber metric, blue distance, pink boundary, green coherence.
# config/metric.yaml
spd_metric:
low_rank:
rank: 4 # Target rank r less or equal D
fit_method: "eigendecomp" # or "svd"
refit_interval: 8 # Iterations between refits
regularization:
delta: 1e-3 # Jitter for numerical stability
adaptive: true # Adjust based on condition number
target_condition: 1e6 # Target cond(M)
whitening:
method: "cholesky" # or "sqrt"
cache_factorization: true # Store W for reuse
validation:
check_spd: true # Verify positive definiteness
warn_threshold: 1e8 # Condition number warning
advanced:
riemannian_updates: false # Use geodesic paths
log_euclidean: false # Work in log-space
parallel_transport: false # Preserve tangent structure{
"session_id": "abc123",
"timestamp": "2025-11-10T14:23:45Z",
"spd_metrics": {
"metric_properties": {
"dimension": 35,
"rank": 4,
"regularization_delta": 0.001,
"condition_number": 1243.56,
"is_well_conditioned": true
},
"distance_computation": {
"mahalanobis_distance": 0.4823,
"euclidean_distance": 0.4163,
"warping_ratio": 1.1585,
"computation_time_ms": 0.23
},
"low_rank_status": {
"current_rank": 4,
"explained_variance": 0.87,
"last_refit_iteration": 16,
"next_refit_iteration": 24
},
"stability_indicators": {
"eigenvalue_min": 0.012,
"eigenvalue_max": 14.92,
"trace": 35.45,
"frobenius_norm": 42.18
}
}
}Compare Euclidean and SPD distances across a grid:
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
import matplotlib.pyplot as plt
def visualize_distance_warping(L, delta=1e-3, grid_size=20):
"""
Visualize how SPD metric warps distances in 2D.
Args:
L: Low-rank factor (2 x r).
delta: Regularization parameter.
grid_size: Grid resolution.
"""
metric = SymbolicManifoldMetric(L, delta)
# Create grid
x = np.linspace(-2, 2, grid_size)
y = np.linspace(-2, 2, grid_size)
X_grid, Y_grid = np.meshgrid(x, y)
# Origin point
origin = np.array([0.0, 0.0])
# Compute distances
d_M = np.zeros((grid_size, grid_size))
d_euclidean = np.zeros((grid_size, grid_size))
for i in range(grid_size):
for j in range(grid_size):
point = np.array([X_grid[i, j], Y_grid[i, j]])
d_M[i, j] = metric.distance(origin, point)
d_euclidean[i, j] = np.linalg.norm(point)
# Plot comparison
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Euclidean distances
im1 = axes[0].contourf(X_grid, Y_grid, d_euclidean, levels=20, cmap='viridis')
axes[0].set_title('Euclidean Distance')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
plt.colorbar(im1, ax=axes[0])
# SPD distances
im2 = axes[1].contourf(X_grid, Y_grid, d_M, levels=20, cmap='plasma')
axes[1].set_title('SPD Warped Distance')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
plt.colorbar(im2, ax=axes[1])
# Ratio
ratio = d_M / (d_euclidean + 1e-10)
im3 = axes[2].contourf(X_grid, Y_grid, ratio, levels=20, cmap='coolwarm')
axes[2].set_title('Warping Ratio (d_M / d_euclidean)')
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
plt.colorbar(im3, ax=axes[2])
plt.tight_layout()
plt.show()
# Example with anisotropic stretching
L_anisotropic = np.array([[2.0, 0.1],
[0.1, 0.5]])
# visualize_distance_warping(L_anisotropic)
print("Visualization code ready (uncomment to display plots)")Apply SPD metrics to educational content analysis:
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
from compitum.metric import SymbolicManifoldMetric
def analyze_essay_clusters(embeddings: np.ndarray, labels: np.ndarray, rank: int = 4):
"""
Cluster essays using an SPD metric.
Args:
embeddings: Essay embeddings (n x D).
labels: Ground truth labels (n,).
rank: Low-rank factor rank.
Returns:
Dictionary with clustering metrics.
"""
# The real constructor is (D, rank, delta) -- it initializes its own
# low-rank factor internally; there is no low_rank_fit() classmethod
# to pre-fit L from data.
D = embeddings.shape[1]
metric = SymbolicManifoldMetric(D=D, rank=rank, delta=1e-3)
# Compute pairwise distances
n = embeddings.shape[0]
distances_M = np.zeros((n, n))
distances_euclidean = np.zeros((n, n))
for i in range(n):
for j in range(i + 1, n):
# distance() returns a (d, sigma) tuple.
d_M, _sigma = metric.distance(embeddings[i], embeddings[j])
d_euc = np.linalg.norm(embeddings[i] - embeddings[j])
distances_M[i, j] = distances_M[j, i] = d_M
distances_euclidean[i, j] = distances_euclidean[j, i] = d_euc
# Analyze within-class vs between-class distances
within_class_M = []
between_class_M = []
within_class_euc = []
between_class_euc = []
for i in range(n):
for j in range(i + 1, n):
if labels[i] == labels[j]:
within_class_M.append(distances_M[i, j])
within_class_euc.append(distances_euclidean[i, j])
else:
between_class_M.append(distances_M[i, j])
between_class_euc.append(distances_euclidean[i, j])
# Compute separation metrics
separation_M = np.mean(between_class_M) / (np.mean(within_class_M) + 1e-10)
separation_euc = np.mean(between_class_euc) / (np.mean(within_class_euc) + 1e-10)
# There is no condition_number() method; compute it directly from the
# real metric_matrix() via eigvalsh if you want to inspect conditioning.
M = metric.metric_matrix()
eigvals = np.linalg.eigvalsh(M)
condition_number = eigvals.max() / eigvals.min()
return {
"separation_spd": separation_M,
"separation_euclidean": separation_euc,
"improvement": separation_M / separation_euc,
"within_mean_spd": np.mean(within_class_M),
"between_mean_spd": np.mean(between_class_M),
"condition_number": condition_number,
}
# Example with synthetic data
np.random.seed(42)
n_samples = 50
D = 20
n_classes = 3
# Create synthetic essay embeddings with class structure
embeddings = []
labels = []
for c in range(n_classes):
center = np.random.randn(D) * 3
class_samples = center + np.random.randn(n_samples // n_classes, D) * 0.5
embeddings.append(class_samples)
labels.extend([c] * (n_samples // n_classes))
embeddings = np.vstack(embeddings)
labels = np.array(labels)
results = analyze_essay_clusters(embeddings, labels, rank=4)
print("Essay Clustering Analysis")
print("=" * 50)
print(f"SPD Separation: {results['separation_spd']:.3f}")
print(f"Euclidean Separation: {results['separation_euclidean']:.3f}")
print(f"Improvement Factor: {results['improvement']:.3f}x")
print(f"Within-class mean (SPD): {results['within_mean_spd']:.3f}")
print(f"Between-class mean (SPD): {results['between_mean_spd']:.3f}")
print(f"Metric condition number: {results['condition_number']:.2e}")Expected Output:
Essay Clustering Analysis
==================================================
SPD Separation: 4.267
Euclidean Separation: 3.845
Improvement Factor: 1.110x
Within-class mean (SPD): 0.892
Between-class mean (SPD): 3.806
Metric condition number: 1.45e+03
| Operation | Time Complexity | Space Complexity | Notes |
|---|---|---|---|
| Build M | O(D times r times r + D squared) | O(D squared) | Low-rank plus regularization |
| Cholesky factorization | O(D cubed) | O(D squared) | One-time per M update |
| Distance computation | O(D squared) | O(D) | Using cached W |
| Low-rank fit | O(n times D squared + D cubed) | O(D squared) | Covariance plus eigen |
| Condition number | O(D cubed) | O(D squared) | Full eigendecomposition |
Dimension Scaling:
| Dimension D | Rank r | Build Time | Distance Time | Recommended Use |
|---|---|---|---|---|
| D less than 20 | D over 2 | Less than 1ms | Less than 0.1ms | Real-time interactive |
| 20 to 50 | 4 to 8 | 1 to 10ms | 0.1 to 0.5ms | Standard applications |
| 50 to 100 | 4 | 10 to 50ms | 0.5 to 2ms | Batch processing |
| D greater than 100 | 4 with PCA | 50 to 200ms | 2 to 10ms | Offline computation |
graph TD
A[Memory Strategy] --> B{Dimension?}
B -->|D less than 50| C[Store full M and W]
B -->|50 less or equal D less than 100| D[Store L and W only]
B -->|D greater or equal 100| E[Compute on demand]
C --> F[Fast access 8D squared bytes]
D --> G[Moderate 4Dr plus 4D squared bytes]
E --> H[Minimal 4Dr bytes]
style A fill:#e3f2fd,stroke:#01579b
style F fill:#e8f5e9,stroke:#1b5e20
style G fill:#fff3e0,stroke:#ef6c00
style H fill:#fff3e0,stroke:#ef6c00
Caption: Memory strategy by dimension. Alt text: Decision tree with 8 nodes showing dimension-based memory strategies, blue start, green fast, amber moderate/minimal.
graph TD
F1[Fundamentals: Read Overview] --> F2[Run Basic Example]
F2 --> F3[Visualize Warping]
F3 --> T1[Tuning: Adjust Rank r]
T1 --> T2[Validate Condition Number]
T2 --> T3[Optimize for Domain]
T3 --> A1[Advanced: Implement Geodesics]
A1 --> A2[Test Parallel Transport]
A2 --> P1[Production: Monitor Stability]
P1 --> P2[Benchmark Performance]
style F1 fill:#e3f2fd,stroke:#01579b
style F2 fill:#e3f2fd,stroke:#01579b
style F3 fill:#e3f2fd,stroke:#01579b
style T1 fill:#e8f5e9,stroke:#1b5e20
style T2 fill:#e8f5e9,stroke:#1b5e20
style T3 fill:#e8f5e9,stroke:#1b5e20
style A1 fill:#fff3e0,stroke:#ef6c00
style A2 fill:#fff3e0,stroke:#ef6c00
style P1 fill:#fce4ec,stroke:#c2185b
style P2 fill:#fce4ec,stroke:#c2185b
Caption: Progressive learning stages. Alt text: Timeline with four sections and 10 steps, blue fundamentals, green tuning, amber advanced, pink production.
Week 1: Fundamentals
- Read documentation: Understand SPD properties and whitening
- Run basic example: Compute d_M for toy embeddings
- Visualize warping: Plot distance contours in 2D
- Compare distances: Compute d_M vs d_euclidean ratios
Week 2: Tuning
- Adjust rank: Test r in {2, 4, 8} for your dimension D
- Validate conditioning: Check cond(M) less than 1e6
- Fit from data: Use low_rank_fit on real embeddings
- Monitor stability: Track eigenvalue spectrum over iterations
Week 3: Advanced
- Implement geodesics: Use geodesic_path for metric interpolation
- Test log-Euclidean: Compare affine-invariant vs log-Euclidean means
- Parallel transport: Preserve tangent vectors along paths
- Production batch: Process 1000 distances, profile performance
Ongoing:
- Document profiles: Contribute optimal (r, delta) pairs to repository
- Report issues: Submit bugs for numerical edge cases
- Expand examples: Add domain-specific applications
- Community engagement: Share findings in discussions
- Configuration: Setting rank, delta, refit intervals in YAML
- Design-Notes: Geometric foundations and motivation
- Output-Format: Distance and metric schemas
- Cookbook: Complete example workflows with visualizations
- Glossary: Definitions for SPD, whitening, Riemannian geometry
Foundational Papers:
- Riemann, B. (1854): "Über die Hypothesen, welche der Geometrie zu Grunde liegen"
- Mahalanobis, P.C. (1936): "On the Generalized Distance in Statistics"
- Amari, S. and Nagaoka, H. (2000): Methods of Information Geometry
Modern Texts:
- Bhatia, R. (2007): Positive Definite Matrices
- Pennec, X., Fillard, P., Ayache, N. (2006): "A Riemannian Framework for Tensor Computing"
- Bonnabel, S. (2013): "Stochastic Gradient Descent on Riemannian Manifolds"
Implementation Guides:
- SciPy Linear Algebra: https://docs.scipy.org/doc/scipy/reference/linalg.html
- PyManOpt: https://pymanopt.org/ (Riemannian optimization)
- Geomstats: https://geomstats.github.io/ (Geometry and statistics on manifolds)
- Repository: https://github.com/PaulTiffany/compitum
- Discussions: https://github.com/PaulTiffany/compitum/discussions
- Issues: https://github.com/PaulTiffany/compitum/issues
- Wiki: https://github.com/PaulTiffany/compitum/wiki
- Examples: https://github.com/PaulTiffany/compitum/tree/main/examples
Q: Why use SPD metrics instead of Euclidean distance?
A: Euclidean distance treats all dimensions equally, which may not reflect the actual structure of your data. SPD metrics learn task-relevant geometry from data—stretching important dimensions, compressing irrelevant ones. This improves clustering, outlier detection, and coherence measures.
Q: How do I choose the rank r?
A: General guidelines:
- Start with r equals 4 for most applications
- For D less than 20: use r equals D over 2
- For D greater than 50: keep r equals 4, use PCA if needed
- Validate by checking explained variance (target 80%+)
Q: What if my M becomes ill-conditioned?
A: Three strategies:
- Increase regularization delta (try 1e-2 instead of 1e-3)
- Reduce rank r to limit flexibility
- Normalize L to have norm less than 10
Always monitor condition number and trigger warnings above 1e8.
Q: When should I refit the metric L?
A: Refit when:
- Embedding distribution shifts significantly
- Every N iterations (default N equals 8)
- Condition number exceeds threshold
- Distance computations become unstable
Use adaptive refitting based on data drift detection.
Q: Can I use SPD metrics with neural network embeddings?
A: Yes! Fit L from a batch of embeddings using low_rank_fit. Update periodically as the network trains. This creates an adaptive distance metric that evolves with your model.
Q: What is whitening and why is it important?
A: Whitening transforms data z equals W times (x minus y) so that the covariance becomes identity. This:
- Decouples correlated dimensions
- Makes kernel density estimation isotropic
- Ensures fair distance comparisons across axes
Always whiten before applying methods that assume isotropy (like KDE).
Last Updated: November 10, 2025 · Compitum v1.0 · Open to Contributions