This approach ("Sine Landmark Reduction") uses linearised trilateration—similar to GPS positioning—against a synthetic "sine skeleton" of landmarks.

The main trade-offs:

It is O(N) and deterministic (solves Ax=b instead of iterative gradient descent).

It forces the topology onto a loop structure, so it is less accurate than UMAP for complex manifolds (like Swiss Rolls), but it guarantees a clean layout for user interfaces.

It can project ~9k points (50 dims) to 3D in about 2 seconds on a laptop CPU. Python implementation and math details are in the post. Happy to answer questions!

Would be nice to see it come back. Would love to browse for books and movies on maps again, rather that getting lists regurgitated at me.

UMAP is not O(n^2) it is O(n log n).

``` import numpy as np import time import matplotlib.pyplot as plt from sklearn.base import BaseEstimator, TransformerMixin from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler, MinMaxScaler from sklearn.datasets import load_digits from sklearn.metrics import pairwise_distances from sklearn.manifold import TSNE

try: import umap HAS_UMAP = True except ImportError: HAS_UMAP = False print("Warning: 'umap-learn' not installed. Comparison will be skipped.")

class SineLandmarkReduction(BaseEstimator, TransformerMixin): def __init__(self, n_components=2, n_landmarks=50, mode='data_derived', # 'sine' or 'data_derived' distance_warping=1.0, random_state=42): self.n_components = n_components self.n_landmarks = n_landmarks self.mode = mode self.p = distance_warping self.random_state = random_state self.rng = np.random.RandomState(random_state)

    def _generate_sine_landmarks(self, n_features):
        """Generates the high-dim 'sine skeleton'."""
        a = self.rng.uniform(0.5, 2.0, n_features)
        omega = self.rng.uniform(0.5, 1.5, n_features)
        phi = self.rng.uniform(0, 2 * np.pi, n_features)
        
        t = np.linspace(0, 2 * np.pi, self.n_landmarks)
        
        L_high = (a[:, None] * np.sin(omega[:, None] * t + phi[:, None])).T
        return L_high

    def fit(self, X, y=None):
        self.scaler = StandardScaler()
        X_scaled = self.scaler.fit_transform(X)
        n_samples, n_features = X_scaled.shape
        
        if self.mode == 'sine':
            self.L_high = self._generate_sine_landmarks(n_features)
            l_min, l_max = self.L_high.min(), self.L_high.max()
            x_min, x_max = X_scaled.min(), X_scaled.max()
            self.L_high = (self.L_high - l_min) / (l_max - l_min) * (x_max - x_min) + x_min
            
        else: # 'data_derived'
            indices = self.rng.choice(n_samples, self.n_landmarks, replace=False)
            self.L_high = X_scaled[indices].copy()

        self.pca_landmarks = PCA(n_components=self.n_components)
        self.L_low = self.pca_landmarks.fit_transform(self.L_high)
        
        self.L0_low = self.L_low[0]
        self.L_others_low = self.L_low[1:]
        
        self.A = 2 * (self.L_others_low - self.L0_low)
        self.A_pinv = np.linalg.pinv(self.A)
        
        self.L0_sq_norm = np.sum(self.L0_low**2)
        self.Li_sq_norms = np.sum(self.L_others_low**2, axis=1)
        
        return self

    def transform(self, X):
        X_scaled = self.scaler.transform(X)
        
        D = pairwise_distances(X_scaled, self.L_high, metric='euclidean')
        
        if self.p != 1.0:
            D = np.power(D, self.p)
            
        D_sq = D**2
        
        d0_sq = D_sq[:, 0:1]
        di_sq = D_sq[:, 1:]
        
        term_dist = d0_sq - di_sq
        term_geom = self.Li_sq_norms - self.L0_sq_norm
        
        B = term_dist + term_geom
        
        Y = np.dot(self.A_pinv, B.T).T
        
        return Y

    print(f"Dataset: Digits (N={X.shape[0]}, D={X.shape[1]})")
    print("-" * 40)

    start = time.time()
    umap_reducer = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=42)
    X_umap = umap_reducer.fit_transform(X)
    umap_time = time.time() - start
    print(f"UMAP Time:  {umap_time:.4f} s")

    start = time.time()
    tsne_reducer = TSNE(n_components=2, perplexity=30, random_state=42)
    X_tsne = tsne_reducer.fit_transform(X)
    tsne_time = time.time() - start
    print(f"t-SNE Time: {tsne_time:.4f} s")

    start = time.time()
    slr = SineLandmarkReduction(n_landmarks=50, mode='data_derived', distance_warping=0.5)
    X_slr = slr.fit_transform(X)
    slr_time = time.time() - start
    print(f"SLR Time:   {slr_time:.4f} s")
    print("-" * 40)
    print(f"SLR vs UMAP Speedup:  {umap_time / slr_time:.1f}x")
    print(f"SLR vs t-SNE Speedup: {tsne_time / slr_time:.1f}x")

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    sc1 = axes[0].scatter(X_umap[:, 0], X_umap[:, 1], c=y, cmap='Spectral', s=5, alpha=0.7)
    axes[0].set_title(f"UMAP\nTime: {umap_time:.3f}s")
    axes[0].axis('off')
    
    sc2 = axes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap='Spectral', s=5, alpha=0.7)
    axes[1].set_title(f"t-SNE\nTime: {tsne_time:.3f}s")
    axes[1].axis('off')
    
    sc3 = axes[2].scatter(X_slr[:, 0], X_slr[:, 1], c=y, cmap='Spectral', s=5, alpha=0.7)
    axes[2].set_title(f"Sine Landmark Reduction (SLR)\nTime: {slr_time:.3f}s")
    axes[2].axis('off')
    
    plt.tight_layout()
    plt.savefig('comparison_plot.png', dpi=150, bbox_inches='tight')
    print("\nPlot saved to comparison_plot.png")
    plt.show()