Repository navigation
Math Kernel Density and Coherence
Quick-Scan Infobox
Aspect Details Target User Educators assessing pattern alignment, developers incorporating trace weights—statistics and machine learning basics assumed. Core Goal Explain how kernel density estimation approximates neighborhood density for coherence, providing a smoothed measure of fit in transformed coordinates. Time to First Insight 3 minutes: Fit on sample traces, compute log density > -2 for alignment. Key Analogy Gathering density: Kernel density estimation averages local points like surveying attendance in areas—clusters show concentration, voids indicate isolation. Compitum Integration Log density in coherence module; adjusts utility, signals edges. Known Limitation High dimensions reduce accuracy due to sparsity; reduce to 3-5 via principal components for reliability.
This structure processes traces and queries through transformation and estimation to yield a coherence value.
graph TD
A[Traces X] --> B[SPD Metric M]
B --> C[Whitening Z equals M inverse sqrt X]
C --> D[KDE Estimator]
D --> E[Smooth Density p-hat z]
F[Query q] --> C
C --> D
D --> G{log p-hat q z}
G --> |greater than -1| H[High: On-Manifold]
G --> |-1 to -3| I[Moderate: Edge]
G --> |less than -3| J[Low: Off-Manifold]
H --> K[Utility Boost]
I --> K
J --> K[Penalty]
style A fill:#fce4ec,stroke:#c2185b
style C fill:#e3f2fd,stroke:#01579b
style D fill:#e8f5e9,stroke:#1b5e20
style H fill:#e8f5e9,stroke:#1b5e20
style I fill:#fff3e0,stroke:#ef6c00
style J fill:#fce4ec,stroke:#c2185b
style K fill:#e8f5e9,stroke:#1b5e20
Caption: Path from inputs to scored utility. Alt text: Graph with 15 nodes showing traces/metrics/whitening/KDE/density/checks/utility flow, pink inputs, blue whitening, green KDE, green high, amber moderate, pink low, green output.
Whitening decouples the metric before estimation. Limitation: Diagram assumes Gaussian kernel; uniform variants alter smoothing.
Coherence evaluates if a query matches a model's local structure: Kernel density estimation computes density p-hat(z) in whitened coordinates z (via SPD inverse M^{-1/2}), giving log p-hat(z) as a balanced indicator of neighborhood consistency.
sequenceDiagram
participant Raw as Traces n x d
participant Metric as SPD M
participant White as Whitening
participant KDE as Density Estimator
participant Query as Query q
participant Score as Coherence Value
Raw->>Metric: Request M
Metric->>White: SPD M dimension d x d
White->>White: Eigendecomposition M equals U Lambda U-transpose
White->>White: Inverse sqrt M inverse sqrt equals U Lambda inverse sqrt U-transpose
Raw->>White: Transform traces
White->>KDE: Z traces equals X M inverse sqrt
Query->>White: Transform query
White->>KDE: z query equals q M inverse sqrt
KDE->>KDE: Fit p-hat z equals sum K z minus z i over h
KDE->>Score: Evaluate log p-hat z query
Score->>Score: Threshold: -1 coherent, -3 incoherent
Score-->>Query: Log evidence
Caption: Sequence from data to score. Alt text: Sequence diagram with six participants and 11 messages, flowing from raw data to final score, green highlighting for fit stage.
Kernel Density Formula:
p-hat(z) = (1 / n h^d) * sum_{i=1 to n} K((z - z_i) / h)
Where:
- K: Kernel (Gaussian standard)
- h: Bandwidth
- n: Points, d: Dimensions
Gaussian Kernel:
K(u) = (2 pi)^{-d/2} * exp(-|u|^2 / 2)
Log Density (Coherence):
log p-hat(z) = -log n - (d/2) log(2 pi h^2) - (1 / 2h^2) min_i |z - z_i|^2
graph TD
A1956[Rosenblatt 1956: Histogram Smoother] --> A1962[Parzen 1962: Kernel Extension]
A1962 --> A1986[Silverman 1986: Bandwidth Optimization]
A1986 --> A1992[Scott 1992: Multivariate Refinements]
A1992 --> A2000[Riemannian KDE 2000s: Compitum Approach]
Caption: Key milestones in kernel density estimation. Alt text: Timeline with five entries from 1956 to 2000s, showing progression of KDE methods.
Innovations:
- Rosenblatt (1956): Histogram continuity foundation
- Parzen (1962): Kernel generalization
- Silverman (1986): Bandwidth optimization via validation
- Scott (1992): Multivariable sparsity handling
Limitation: Timeline focuses on KDE fundamentals; Riemannian extensions developed post-2000.
Mean Squared Error:
- Bias: O(h^2) (decreases with h)
- Variance: O(1/(n h^d)) (increases with h)
- Optimal: Balance determined by validation
Compitum Application:
- High coherence (>-1 nat) correlates with low variance (0.05)
- Threshold -3 triggers utility penalties
- Whitening uniformizes h across geometries
Limitation: Assumes smooth underlying density; discontinuous data requires robust kernels.
graph LR
subgraph Input[Input Space]
A[Anisotropic SPD Metric M]
B[Distorted Coordinates]
end
subgraph Transform[Whitening]
C[Eigendecomposition: M equals U Lambda U-transpose]
D[Inverse Square Root: M inverse sqrt equals U Lambda inverse sqrt U-transpose]
E[Transform: Z equals X M inverse sqrt]
end
subgraph Output[Whitened Space]
F[Isotropic Identity Metric]
G[Even Distances]
end
A --> C
B --> E
C --> D
D --> E
E --> F
E --> G
style Input fill:#fce4ec,stroke:#c2185b
style Transform fill:#e3f2fd,stroke:#01579b
style Output fill:#e8f5e9,stroke:#1b5e20
Caption: Transformation from metric to isotropic space. Alt text: Graph with three subgraphs and 9 nodes, showing pink input, blue transform, and green output stages.
graph TD
START[Select Bandwidth h] --> Q1{Validation Set?}
Q1 -->|Yes| CV[Cross-Validation]
CV --> GRID[h in 0.05-0.5 range]
GRID --> EV[Maximize Log-Likelihood]
EV --> DONE[Optimal h]
Q1 -->|No| RULE{Known n, sigma?}
RULE -->|Yes| SIL[Silverman Rule]
SIL --> FORM[h equals 1.06 sigma n to power -0.2]
FORM --> DONE
RULE -->|No| DEF[h equals 0.1 default]
DEF --> TEST{Visual Check?}
TEST -->|Spiky| INC[Increase h by 50 percent]
TEST -->|Flat| DEC[Decrease h by 50 percent]
TEST -->|Fine| DONE
INC --> DONE
DEC --> DONE
style START fill:#e3f2fd,stroke:#01579b
style DONE fill:#e8f5e9,stroke:#1b5e20
style CV fill:#fff3e0,stroke:#ef6c00
style SIL fill:#e8f5e9,stroke:#1b5e20
style DEF fill:#fff3e0,stroke:#ef6c00
Caption: Decision tree for bandwidth selection. Alt text: Graph with start node branching to validation/rule/default paths, 18 nodes total, blue start, green completion, amber CV, green Silverman, amber default.
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
class CoherenceEstimator:
"""
Density-based coherence estimation in whitened space.
Estimates local density to assess query alignment with trace manifold.
"""
def __init__(self, bandwidth: float = 0.1, threshold_coherent: float = -1.0,
threshold_incoherent: float = -3.0):
"""
Initialize coherence estimator.
Args:
bandwidth: Smoothing parameter h (uses Silverman if None).
threshold_coherent: Above this value indicates alignment.
threshold_incoherent: Below this value indicates misalignment.
"""
self.h = bandwidth
self.threshold_high = threshold_coherent
self.threshold_low = threshold_incoherent
self.kde = None
self.M_inv_sqrt = None
def _whiten(self, X: np.ndarray, M: np.ndarray) -> np.ndarray:
"""
Apply isotropic transformation via metric inverse square root.
Args:
X: Points (n x d).
M: SPD matrix (d x d).
Returns:
Z: Whitened coordinates (n x d).
"""
eigenvalues, U = np.linalg.eigh(M)
s_inv = np.where(eigenvalues > 1e-10, 1/np.sqrt(eigenvalues), 0)
self.M_inv_sqrt = U @ np.diag(s_inv) @ U.T
return X @ self.M_inv_sqrt.T
def fit(self, traces: np.ndarray, M: np.ndarray) -> 'CoherenceEstimator':
"""
Fit KDE on whitened traces.
Args:
traces: Training points (n x d).
M: Metric matrix (d x d).
Returns:
Self for chaining.
"""
Z = self._whiten(traces, M)
self.kde = gaussian_kde(Z.T, bw_method=self.h)
return self
def score(self, query: np.ndarray) -> tuple[float, str]:
"""
Compute log density and coherence label.
Args:
query: Query point (d,) or (1 x d).
Returns:
Tuple of (log probability, label).
"""
if self.kde is None:
raise ValueError("Must call fit() before scoring.")
z_query = query.reshape(1, -1) @ self.M_inv_sqrt.T
log_p = np.log(self.kde(z_query.T)[0] + 1e-10)
if log_p > self.threshold_high:
label = "high"
elif log_p > self.threshold_low:
label = "moderate"
else:
label = "low"
return log_p, label
def visualize_density(self, traces: np.ndarray, query: np.ndarray,
M: np.ndarray, grid_size: int = 50):
"""
Generate 2D density heatmap (requires d=2).
Args:
traces: Training points (n x 2).
query: Query point (2-vector).
M: Metric matrix (2 x 2).
grid_size: Grid resolution for visualization.
"""
if traces.shape[1] != 2:
raise ValueError("Visualization requires 2D data.")
Z_traces = self._whiten(traces, M)
z_query = self._whiten(query.reshape(1, -1), M)
x_min, x_max = Z_traces[:, 0].min() - 1, Z_traces[:, 0].max() + 1
y_min, y_max = Z_traces[:, 1].min() - 1, Z_traces[:, 1].max() + 1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, grid_size),
np.linspace(y_min, y_max, grid_size))
grid_points = np.c_[xx.ravel(), yy.ravel()]
densities = np.exp(self.kde.logpdf(grid_points.T))
densities = densities.reshape(xx.shape)
plt.figure(figsize=(10, 8))
plt.contourf(xx, yy, densities, levels=20, cmap='viridis', alpha=0.8)
plt.colorbar(label='Density p-hat(z)')
plt.scatter(Z_traces[:, 0], Z_traces[:, 1], c='white',
s=30, edgecolors='black', label='Traces')
plt.scatter(z_query[:, 0], z_query[:, 1], c='red',
s=200, marker='*', edgecolors='black', label='Query')
log_p_query = self.kde.logpdf(z_query.T)[0]
plt.contour(xx, yy, np.log(densities + 1e-10),
levels=[self.threshold_low, self.threshold_high, log_p_query],
colors=['red', 'orange', 'yellow'], linewidths=2)
plt.xlabel('Z_1')
plt.ylabel('Z_2')
plt.title(f'Density Map (h={self.h:.3f})')
plt.legend()
plt.grid(alpha=0.3)
plt.show()
# Example Usage
np.random.seed(42)
traces = np.random.multivariate_normal([0, 0], [[1, 0.3], [0.3, 1]], 100)
M = np.eye(2)
q_coherent = np.array([0.1, 0.1])
q_incoherent = np.array([3.5, 3.5])
estimator = CoherenceEstimator(bandwidth=0.2)
estimator.fit(traces, M)
log_p_coh, label_coh = estimator.score(q_coherent)
log_p_inc, label_inc = estimator.score(q_incoherent)
print("Coherence Analysis")
print("=" * 50)
print(f"Coherent query: log p-hat = {log_p_coh:.3f} ({label_coh})")
print(f"Incoherent query: log p-hat = {log_p_inc:.3f} ({label_inc})")
print("\nInterpretation:")
print(f" > -1.0: High (on-manifold)")
print(f" -1 to -3: Moderate (edge)")
print(f" < -3.0: Low (off-manifold)")
# Uncomment to visualize:
# estimator.visualize_density(traces, q_coherent, M)Expected Output:
Coherence Analysis
==================================================
Coherent query: log p-hat = -0.847 (high)
Incoherent query: log p-hat = -5.231 (low)
Interpretation:
> -1.0: High (on-manifold)
-1 to -3: Moderate (edge)
< -3.0: Low (off-manifold)
Note: Values depend on random seed; visualization shows density contours if uncommented. Limitation: Visualization limited to 2D; generalize via dimension reduction or slicing.
graph TD
subgraph Space[Raw Space]
A[Dense Cluster]
B[Sparse Areas]
C[Central Group]
end
subgraph Estimation[Smoothing Process]
D[Local Counts]
E[Kernel Average]
F[Density Map]
end
subgraph Query[Query Placement]
G[In Cluster: High]
H[Near Edge: Medium]
I[Outside: Low]
end
A --> D
B --> D
C --> D
D --> E
E --> F
F --> G
F --> H
F --> I
style Space fill:#fce4ec,stroke:#c2185b
style Estimation fill:#e3f2fd,stroke:#01579b
style G fill:#e8f5e9,stroke:#1b5e20
style H fill:#fff3e0,stroke:#ef6c00
style I fill:#fce4ec,stroke:#c2185b
Caption: Transformation from raw data to query placement. Alt text: Graph with three subgraphs and 12 nodes, showing pink space, blue estimation, green high, amber medium, pink low placement.
Essay Alignment:
Embed essays into vector space; fit KDE on exemplar traces. New essays score by density proximity.
Use: "High density indicates strong thematic match; low density suggests redirection to core ideas."
Outlier Identification:
Low density values flag deviations from established patterns.
Use: "Group outliers for targeted guidance and support."
Adaptive Grouping:
Cluster students by density regions for differentiated instruction.
graph LR
subgraph Tight[Narrow Bandwidth]
T0[h equals 0.01]
A1[Fragmented Detail]
A2[Noise Amplification]
A3[Uneven Estimates]
end
subgraph Balanced[Appropriate Bandwidth]
B0[h equals 0.1]
B1[Clear Forms]
B2[Balanced Detail]
B3[Consistent Scores]
end
subgraph Wide[Broad Bandwidth]
C0[h equals 1.0]
C1[Over-Smoothed]
C2[Lost Edges]
C3[Uniform Levels]
end
Tight -.-> Balanced
Balanced -.-> Wide
style Tight fill:#fce4ec,stroke:#c2185b
style Balanced fill:#e8f5e9,stroke:#1b5e20
style Wide fill:#e3f2fd,stroke:#01579b
Caption: Effects of bandwidth selection. Alt text: Graph with three subgraphs and 9 nodes, showing pink tight, green balanced, blue wide bandwidth effects.
Montessori Parallel: Like color tablets revealing hue clusters, density estimation shows natural groupings—a child matches by local proximity (kernel averaging), adjusting scale (bandwidth) for accurate pairing, learning discernment through neighborhood relationships.
Limitation: Parallels simplify concepts; validate with domain-specific data.
graph LR
subgraph Inputs[Input Components]
I1[Traces X: n x d]
I2[Metric M: d x d]
I3[Query q: d dimensions]
end
subgraph Processing[Core Processing]
P1[Whitening Transform]
P2[KDE Fitting]
P3[Density Evaluation]
end
subgraph Outputs[Output Signals]
O1[Log Evidence]
O2[Coherence Label]
O3[Utility Weight]
end
I1 --> P1
I2 --> P1
I3 --> P1
P1 --> P2
P2 --> P3
P3 --> O1
O1 --> O2
O2 --> O3
style Inputs fill:#e3f2fd,stroke:#01579b
style Processing fill:#e8f5e9,stroke:#1b5e20
style Outputs fill:#fff3e0,stroke:#ef6c00
Caption: Flow from inputs to weighted outputs. Alt text: Graph with three subgraphs and 15 nodes, showing blue inputs, green processing, amber outputs.
| Element | Definition | Key Choice | Range | Interpretation | Compitum Field |
|---|---|---|---|---|---|
| Kernel K | Influence shape | Gaussian (default) | N/A | Spread of averaging | Fixed in implementation |
| Bandwidth h | Smoothing scale | Validation or Silverman | 0.05–0.5 | Bias-variance balance | coherence.bandwidth |
| Whitening | z = M^{-1/2} x | Eigendecomposition | Always applied | Metric uniformization | Preprocessing step |
| Log Evidence | log p-hat(z) | Threshold-based | -infinity to 0 | Alignment measure | utility.coherence_score |
| High Threshold | Coherent limit | Domain-tuned | -1.0 (typical) | On-structure indicator | config.threshold_high |
| Low Threshold | Incoherent limit | Domain-tuned | -3.0 (typical) | Off-structure indicator | config.threshold_low |
graph LR
A[Log Evidence] --> B{Score Range?}
B -->|greater than -0.5| C[Exceptional: Core Structure]
B -->|-0.5 to -1.0| D[High: Aligned]
B -->|-1.0 to -2.0| E[Moderate: Perimeter]
B -->|-2.0 to -3.0| F[Low: Dubious]
B -->|less than -3.0| G[Very Low: Divergent]
C --> H[Weight multiplier 2.0]
D --> I[Weight multiplier 1.0]
E --> J[Weight multiplier 0.5]
F --> K[Weight multiplier 0.2]
G --> L[Weight multiplier 0.0]
style C fill:#e8f5e9,stroke:#1b5e20
style D fill:#e8f5e9,stroke:#1b5e20
style E fill:#fff3e0,stroke:#ef6c00
style F fill:#fff3e0,stroke:#ef6c00
style G fill:#fce4ec,stroke:#c2185b
Caption: Score ranges and corresponding utility weights. Alt text: Graph showing log evidence branching to ranges and weights, 15 nodes total, green high, amber moderate, pink low categories.
graph TD
START[Problem Detected] --> Q1{Density Quality Issue?}
Q1 -->|Flat density| Q2{h greater than 0.5?}
Q2 -->|Yes| FIX1[Reduce h to 0.1; ensure n greater than 20]
Q2 -->|No| FIX2[Add more traces]
Q1 -->|Spiky density| Q3{h less than 0.05?}
Q3 -->|Yes| FIX3[Increase h to 0.2; clip p min at 1e-10]
Q3 -->|No| FIX4[Filter outliers]
Q1 -->|Skewed results| Q4{Whitening Applied?}
Q4 -->|No| FIX5[Apply M inverse sqrt transform]
Q4 -->|Yes| Q5{M Numerically Stable?}
Q5 -->|No| FIX6[Regularize: M plus lambda I]
Q5 -->|Yes| FIX7[Verify metric construction]
Q1 -->|NaN values| Q6{Log of Zero?}
Q6 -->|Yes| FIX8[Add epsilon: p min 1e-10]
Q6 -->|No| FIX9[Check eigenvalues greater than 0]
style START fill:#fce4ec,stroke:#c2185b
style FIX1 fill:#e8f5e9,stroke:#1b5e20
style FIX2 fill:#e8f5e9,stroke:#1b5e20
style FIX3 fill:#e8f5e9,stroke:#1b5e20
style FIX4 fill:#e8f5e9,stroke:#1b5e20
style FIX5 fill:#e8f5e9,stroke:#1b5e20
style FIX6 fill:#e3f2fd,stroke:#01579b
style FIX7 fill:#e8f5e9,stroke:#1b5e20
style FIX8 fill:#e8f5e9,stroke:#1b5e20
style FIX9 fill:#fff3e0,stroke:#ef6c00
Caption: Resolution path for common issues. Alt text: Graph showing problem root branching to density quality checks, 20 nodes total, pink start, green fixes, blue regularization, amber numerical fixes.
| Issue | Symptom | Cause | Fix | Parallel |
|---|---|---|---|---|
| Over-Smooth | Flat p-hat near 0 | h > 0.5 or n small | Set h=0.05; ensure n>50 | Over-blurred photograph |
| Under-Smooth | Spikes, NaN in log | h < 0.02 | Set h=0.2; clip p >= 1e-10 | Grainy, pixelated scan |
| Anisotropic Bias | Distorted densities | No whitening applied | Apply M^{-1/2} transform | Viewing through skewed lens |
| Dimension Curse | High variance in d>5 | Data sparsity | Reduce to 3D via PCA; use product kernels | Vast empty warehouse |
| Singular M | Eigendecomp crash | Ill-conditioned metric | Regularize: M + lambda*I (lambda=1e-6) | Broken measurement tool |
| Threshold Mismatch | All high or all low | Domain mismatch | Tune on validation set | Ill-calibrated gauge |
Limitation: Fixes address standard cases; log unique situations for debugging.
graph TD
START[Begin Tuning] --> SPLIT[Split: Train 80 percent, Val 20 percent]
SPLIT --> GRID[Grid: h from 0.05 to 0.5, step 0.05]
GRID --> LOOP[For each h value]
LOOP --> FIT[Fit KDE on Train]
FIT --> EVAL[Evaluate Log-Likelihood on Val]
EVAL --> STORE[Store h and score]
STORE --> MORE{More h values?}
MORE -->|Yes| LOOP
MORE -->|No| MAX[Select h with Max Score]
MAX --> VIS[Visualize Density]
VIS --> CHECK{Qualitatively Good?}
CHECK -->|Yes| DONE[Deploy h]
CHECK -->|No| REVISE[Adjust Grid Range]
REVISE --> GRID
style START fill:#e3f2fd,stroke:#01579b
style GRID fill:#fff3e0,stroke:#ef6c00
style MAX fill:#e8f5e9,stroke:#1b5e20
style DONE fill:#e8f5e9,stroke:#1b5e20
style REVISE fill:#fff3e0,stroke:#ef6c00
Caption: Systematic bandwidth selection process. Alt text: Flowchart with 18 nodes showing start to split/grid/evaluation/max/visualization/check cycle, blue start, amber grid, green max/done, amber revision.
Profiles:
| Use Case | h | High Threshold | Low Threshold | Notes |
|---|---|---|---|---|
| Dense Data | 0.1 | -0.5 | -2.0 | n > 100 samples |
| Standard | 0.2 | -1.0 | -3.0 | Default configuration |
| Sparse Data | 0.3 | -1.5 | -4.0 | n = 20-50 samples |
| High Dimension | 0.15 | -2.0 | -5.0 | After PCA to d=3-5 |
| Initial Survey | 0.25 | -1.2 | -3.5 | Exploratory phase |
Silverman Bandwidth Calculator:
Open on GitHub | Open in nbviewer | Launch in Binder
def silverman_bandwidth(X: np.ndarray) -> float:
"""
Compute optimal bandwidth via Silverman's rule of thumb.
Args:
X: Data matrix (n x d).
Returns:
Optimal bandwidth h.
"""
n, d = X.shape
sigma = np.std(X, axis=0).mean()
h = 1.06 * sigma * n ** (-1 / (d + 4))
return h
# Example
traces = np.random.randn(100, 3)
h_opt = silverman_bandwidth(traces)
print(f"Optimal h = {h_opt:.3f}")Expected Output:
Optimal h = 0.123
Limitation: Rule assumes Gaussian distribution; Epanechnikov kernel requires 0.9x adjustment.
graph TD
A[Start Checklist] --> B{Visual Quality?}
B -->|Pass| C{Numeric Stability?}
B -->|Fail| D[Adjust h by 50 percent]
C -->|Pass| E{Threshold Alignment?}
C -->|Fail| F[Apply clipping or whitening]
E -->|Pass| G{Production Ready?}
E -->|Fail| H[Retune thresholds]
G -->|Pass| I[Deploy System]
G -->|Fail| J[Collect more traces]
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 checklist workflow. Alt text: Graph showing checklist flow through visual/numeric/threshold/production branches, 15 nodes, blue start, green deploy, amber revision steps.
Checklist Items:
| Phase | Test | Pass Criteria | Fail Action |
|---|---|---|---|
| Visual | Smooth heatmap | No excessive peaks or valleys | Adjust h by 50% |
| Visual | Labels accurate | Manual review matches | Retune parameters |
| Numeric | No NaN or infinity | All values finite | Clip p >= 1e-10 |
| Numeric | Eigenvalues positive | min(eigenvalues) > 1e-8 | Regularize M |
| Threshold | 60-80% moderate | Balanced distribution | Adjust thresholds |
| Threshold | Less than 5% errors | Versus known labels | Widen or narrow range |
| Production | Aligned queries routed well | Utility gain observed | Refit KDE |
| Production | Misaligned queries flagged | Boundary detection works | Lower low threshold |
graph TD
subgraph HighD[High Dimensionality d greater than 10]
A[Traces X: n x 50]
B[Sparsity Curse]
end
subgraph PCA[Projection Stage]
C[Compute Covariance]
D[Extract Top k equals 3 Components]
E[Reduced X prime]
end
subgraph KDE[Low Dimension KDE]
F[Fit on Reduced Data]
G[Reliable p-hat Estimates]
H[Valid Coherence Score]
end
A --> C
B -.-> C
C --> D
D --> E
E --> F
F --> G
G --> H
style A fill:#fce4ec,stroke:#c2185b
style B fill:#fce4ec,stroke:#c2185b
style E fill:#e3f2fd,stroke:#01579b
style H fill:#e8f5e9,stroke:#1b5e20
Caption: Dimension reduction before density estimation. Alt text: Graph with three subgraphs and 10 nodes, showing high-d problem to PCA to KDE stages, pink high-d, blue reduced, green valid output.
Why PCA First:
- d=10 volume: 2^10 = 1024; n=100 gives density ≈ 0.098
- d=3 volume: 2^3 = 8; density ≈ 12.5 (128x improvement)
Open on GitHub | Open in nbviewer | Launch in Binder
from sklearn.decomposition import PCA
def fit_kde_with_pca(traces: np.ndarray, M: np.ndarray,
n_components: int = 3) -> tuple:
"""
Fit KDE after PCA dimensionality reduction.
Args:
traces: High-dimensional data (n x d).
M: Metric matrix (d x d).
n_components: Target dimensionality.
Returns:
Tuple of (estimator, pca_transformer).
"""
estimator = CoherenceEstimator()
Z = estimator._whiten(traces, M)
pca = PCA(n_components=n_components)
Z_reduced = pca.fit_transform(Z)
estimator.kde = gaussian_kde(Z_reduced.T)
estimator.pca = pca
variance_explained = pca.explained_variance_ratio_.sum()
print(f"Variance retained: {variance_explained:.2%}")
return estimator, pca
# Example
traces_high_dim = np.random.randn(100, 20)
M_high_dim = np.eye(20)
est, pca = fit_kde_with_pca(traces_high_dim, M_high_dim, n_components=5)Expected Output:
Variance retained: 75%
Limitation: PCA loses some variance; recommended only for d > 10.
graph LR
subgraph Kernels[Kernel Functions]
G[Gaussian: Smooth, Infinite Support]
E[Epanechnikov: Efficient, Compact]
T[Triangular: Simple, Linear]
U[Uniform: Box, Discontinuous]
end
subgraph Properties[Trade-offs]
P1[Smoothness]
P2[Computation Speed]
P3[Edge Effects]
end
G -.->|Highest| P1
G -.->|Slower| P2
G -.->|Better| P3
E -.->|Medium| P1
E -.->|Fastest| P2
E -.->|Medium| P3
T -.->|Lower| P1
T -.->|Fast| P2
T -.->|Worse| P3
U -.->|None| P1
U -.->|Very Fast| P2
U -.->|Worst| P3
style G fill:#e8f5e9,stroke:#1b5e20
style E fill:#e8f5e9,stroke:#1b5e20
style T fill:#fff3e0,stroke:#ef6c00
style U fill:#fce4ec,stroke:#c2185b
Caption: Kernel function trade-offs. Alt text: Graph with kernels subgraph and properties subgraph, 15 nodes total, green Gaussian/Epanechnikov, amber Triangular, pink Uniform.
Selection Guide:
- Gaussian: Smooth transitions, no discontinuities, infinite tails
- Epanechnikov: Optimal MSE, compact support, 30% faster computation
- Triangular: Basic option, good for initial testing
- Uniform: Avoid (introduces jumps and discontinuities)
Limitation: Gaussian is default; switch to Epanechnikov for performance optimization.
graph TD
F1[Fundamentals: Read Documentation] --> F2[Run Sample Code]
F2 --> F3[Visualize Density Maps]
F3 --> T1[Tuning: Experiment with h]
T1 --> T2[Validate Threshold Limits]
T2 --> T3[Optimize for Your Data]
T3 --> A1[Advanced: Add PCA for High-d]
A1 --> A2[Test Alternative Kernels]
A2 --> P1[Production: Track Coherence Metrics]
P1 --> P2[Benchmark Utility Improvements]
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 from fundamentals to production. Alt text: Timeline with four sections and 12 steps, blue fundamentals, green tuning, amber advanced, pink production.
Action Items:
Week 1: Fundamentals
- Fit KDE: Use sample traces from Cookbook with h in {0.05, 0.1, 0.2}, plot log p(z)
- Visualize: Create 2D heatmap with aligned and misaligned queries
- Whitening test: Compare results with and without M^{-1/2} on skewed data
Week 2: Tuning
- Adjust coherence weight: Set beta_s=0.15 in config, route queries, note utility shifts
- Threshold sweep: Test [-0.5, -1.0, -1.5] for high threshold, observe utility changes
- Cross-validation: Hold out 20% of data, maximize cross-log-likelihood
Week 3: Advanced
- High-dimension test: Apply PCA on d>10 data, compare KDE performance
- Kernel comparison: Implement Epanechnikov in scipy, benchmark speed
- Production batch: Process 100 samples, analyze coherence score distribution
Ongoing:
- Document profiles: Contribute optimal (h, threshold) pairs for your domain to repository
- Report issues: Submit bugs for singular M handling, NaN edge cases
- Expand documentation: Add examples for text embeddings, images, time series
Expand for Full Derivation
Step 1: Whitening Transformation
M = U Lambda U^T (eigendecomposition)
M^{-1/2} = U Lambda^{-1/2} U^T
Z = X M^{-1/2}
Property: |z_i - z_j|_Identity = |x_i - x_j|_M
Step 2: Kernel Density Estimation
p-hat(z) = (1 / n h^d) * sum_{i=1 to n} K((z - z_i) / h)
For Gaussian kernel:
K(u) = (2 pi)^{-d/2} * exp(-|u|^2 / 2)
p-hat(z) = (1 / (n (2 pi h^2)^{d/2})) * sum_{i=1 to n} exp(-|z - z_i|^2 / (2 h^2))
Step 3: Log Evidence
log p-hat(z) = -log n - (d/2) log(2 pi h^2) + log(sum_{i=1 to n} exp(-|z - z_i|^2 / (2 h^2)))
Nearest neighbor approximation:
log p-hat(z) ≈ -log n - (d/2) log(2 pi h^2) - r_min^2 / (2 h^2)
where r_min = min_i |z - z_i|
Step 4: Asymptotic Analysis
- Bias: O(h^2) * integral(u^2 K(u) du) / 2
- Variance: O(1/(n h^d)) * p(z) * integral(K^2(u) du)
- MSE: O(h^4) + O(1/(n h^d))
Optimal bandwidth scales as h ~ n^{-1/(d+4)}
For Gaussian in d=2:
h = (4 / ((d+2) n))^{1/(d+4)} * sigma ≈ 1.06 * sigma * n^{-1/6}
Step 5: Threshold Interpretation
By empirical quantiles:
| Quantile | Log Evidence | Label | Meaning |
|---|---|---|---|
| Top 20% | > -0.5 | Exceptional | Central manifold structure |
| 20-60% | -0.5 to -1.5 | High | Well-aligned |
| 60-85% | -1.5 to -3.0 | Moderate | Perimeter region |
| 85-95% | -3.0 to -5.0 | Low | Divergent |
| Bottom 5% | < -5.0 | Very Low | Anomaly |
Utility Integration:
U_coherence = beta_s * min(max(log p-hat(z), -5), 0)
where beta_s in [0.1, 0.3] contributes to free energy calculation.
- Configuration: Setting h and thresholds in YAML
- Design-Notes: Utility and coherence integration details
- Output-Format: Evidence value schema
- Cookbook: Complete example loops with visualizations
- Glossary: Definitions for whitened space, SPD, log evidence
Foundational Papers:
- Rosenblatt (1956): "Remarks on Some Nonparametric Estimates of a Density Function"
- Parzen (1962): "On Estimation of a Probability Density Function and Mode"
- Silverman (1986): Density Estimation for Statistics and Data Analysis
- Scott (1992): Multivariate Density Estimation
Advanced Topics:
- Wand and Jones (1995): Kernel Smoothing
- Gramacki (2018): Nonparametric Kernel Density Estimation and Its Computational Aspects
- Pelletier (2005): "Kernel Density Estimation on Riemannian Manifolds"
Implementation Guides:
- SciPy gaussian_kde: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html
- Scikit-learn KernelDensity: https://scikit-learn.org/stable/modules/density.html
- 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
Last Updated: November 10, 2025 · Compitum v1.0 · Open to Contributions