Isomap: Isometric Feature Mapping

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

import numpy as np
from scipy.spatial.distance import cdist
 
def random_landmarks(X, n_landmarks, seed=None):
    """
    Select landmarks uniformly at random.
    
    Parameters:
        X: (N, D) array of data points
        n_landmarks: Number of landmarks L
        seed: Random seed for reproducibility
    
    Returns:
        landmark_indices: Array of L indices into X
    """
    rng = np.random.RandomState(seed)
    return rng.choice(X.shape[0], size=n_landmarks, replace=False)
 
def maxmin_landmarks(X, n_landmarks, seed=None):
    """
    Select landmarks using MaxMin (Gonzalez) algorithm.
    
    Iteratively selects the point farthest from all current landmarks.
    Provides good coverage of the data distribution.
    
    Parameters:
        X: (N, D) array of data points
        n_landmarks: Number of landmarks L
        seed: Random seed for first landmark
    
    Returns:
        landmark_indices: Array of L indices into X
    """
    N = X.shape[0]
    rng = np.random.RandomState(seed)
    
    # Initialize with random first landmark
    landmarks = [rng.randint(N)]
    
    # Track minimum distance to any landmark for each point
    min_distances = np.full(N, np.inf)
    
    for _ in range(n_landmarks - 1):
        # Update min distances with last added landmark
        last_landmark = landmarks[-1]
        dists = np.linalg.norm(X - X[last_landmark], axis=1)
        min_distances = np.minimum(min_distances, dists)
        
        # Exclude existing landmarks
        min_distances[landmarks] = -np.inf
        
        # Select point with maximum min-distance
        next_landmark = np.argmax(min_distances)
        landmarks.append(next_landmark)
    
    return np.array(landmarks)
 
def kmeansplusplus_landmarks(X, n_landmarks, seed=None):
    """
    Select landmarks using k-means++ initialization.
    
    Probabilistically selects points proportional to squared distance
    from nearest existing landmark.
    
    Parameters:
        X: (N, D) array of data points
        n_landmarks: Number of landmarks L
        seed: Random seed
    
    Returns:
        landmark_indices: Array of L indices into X
    """
    N = X.shape[0]
    rng = np.random.RandomState(seed)
    
    # Initialize with random first landmark
    landmarks = [rng.randint(N)]
    
    for _ in range(n_landmarks - 1):
        # Compute squared distances to nearest landmark
        landmark_points = X[landmarks]
        dists = cdist(X, landmark_points, metric='sqeuclidean')
        min_sq_dists = dists.min(axis=1)
        
        # Zero out existing landmarks
        min_sq_dists[landmarks] = 0
        
        # Sample proportional to squared distance
        probs = min_sq_dists / min_sq_dists.sum()
        next_landmark = rng.choice(N, p=probs)
        landmarks.append(next_landmark)
    
    return np.array(landmarks)
 
def compare_landmark_strategies(X, n_landmarks, n_trials=10):
    """
    Compare landmark selection strategies by measuring coverage.
    
    Coverage is measured as the average distance from each point
    to its nearest landmark (lower is better).
    """
    N = X.shape[0]
    results = {}
    
    for name, func in [('Random', random_landmarks),
                       ('MaxMin', maxmin_landmarks),
                       ('K-means++', kmeansplusplus_landmarks)]:
        coverages = []
        for trial in range(n_trials):
            landmarks = func(X, n_landmarks, seed=trial)
            landmark_points = X[landmarks]
            dists = cdist(X, landmark_points)
            min_dists = dists.min(axis=1)
            coverages.append(min_dists.mean())
        
        results[name] = {
            'mean_coverage': np.mean(coverages),
            'std_coverage': np.std(coverages)
        }
    
    return results

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

import numpy as np
from scipy.linalg import eigh
 
def landmark_mds(D_landmarks, D_all, n_components):
    """
    Perform Landmark MDS embedding.
    
    Parameters:
        D_landmarks: (L, L) geodesic distance matrix between landmarks
        D_all: (L, N) geodesic distances from landmarks to all points
        n_components: Target embedding dimension d
    
    Returns:
        embedding: (N, n_components) embedded coordinates
        eigenvalues: Eigenvalues (for explained variance)
    """
    L, N = D_all.shape
    
    # Phase 1: Embed landmarks using classical MDS
    # Step 1: Square the distances
    D_sq = D_landmarks ** 2
    
    # Step 2: Double centering
    # B = -0.5 * H @ D_sq @ H, where H = I - (1/L) * 1 @ 1.T
    row_mean = D_sq.mean(axis=1, keepdims=True)
    col_mean = D_sq.mean(axis=0, keepdims=True)
    total_mean = D_sq.mean()
    B = -0.5 * (D_sq - row_mean - col_mean + total_mean)
    
    # Step 3: Eigendecomposition (get top eigenvalues/vectors)
    # Use 'eigh' for symmetric matrices; returns eigenvalues in ascending order
    eigenvalues, eigenvectors = eigh(B)
    
    # Select top n_components (largest eigenvalues)
    idx = np.argsort(eigenvalues)[::-1][:n_components]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    
    # Handle potential negative eigenvalues (due to non-Euclidean distances)
    eigenvalues = np.maximum(eigenvalues, 0)
    
    # Step 4: Landmark coordinates
    Y_landmarks = eigenvectors * np.sqrt(eigenvalues)
    
    # Phase 2: Embed all points using Nyström extension
    # Compute squared distances from all points to landmarks
    D_all_sq = D_all ** 2  # (L, N)
    
    # Column mean of landmark squared distance matrix
    d_bar = D_sq.mean(axis=1)  # (L,)
    
    # Nyström extension: y_n = (1/2) * Lambda^(-1/2) @ V.T @ (d_bar - d_n^2)
    # Rearranged: embedding = (D_all_sq - d_bar[:, None]).T @ V @ Lambda^(-1/2) * (-0.5)
    
    # For numerical stability, use pseudo-inverse
    Lambda_inv_sqrt = np.where(eigenvalues > 1e-10, 
                                1.0 / np.sqrt(eigenvalues), 
                                0)
    
    # embedding = -0.5 * (D_all_sq.T - d_bar) @ eigenvectors @ np.diag(Lambda_inv_sqrt)
    embedding = -0.5 * (D_all_sq.T - d_bar) @ (eigenvectors * Lambda_inv_sqrt)
    
    return embedding, eigenvalues
 
def explained_variance_ratio(eigenvalues):
    """Compute explained variance ratio for each component."""
    eigenvalues = np.maximum(eigenvalues, 0)  # Handle negative eigenvalues
    total = eigenvalues.sum()
    if total > 0:
        return eigenvalues / total
    return np.zeros_like(eigenvalues)
 
class LandmarkIsomap:
    """
    Complete Landmark Isomap implementation.
    
    Combines:
    1. k-NN graph construction
    2. Dijkstra shortest paths from landmarks
    3. Landmark MDS embedding
    """
    
    def __init__(self, n_neighbors=5, n_landmarks=100, n_components=2, 
                 landmark_method='maxmin'):
        self.n_neighbors = n_neighbors
        self.n_landmarks = n_landmarks
        self.n_components = n_components
        self.landmark_method = landmark_method
        
    def fit_transform(self, X):
        """Fit and transform data."""
        from scipy.sparse.csgraph import dijkstra
        
        N = X.shape[0]
        L = min(self.n_landmarks, N)
        
        # Step 1: Build k-NN graph
        adjacency = build_knn_graph_sklearn(X, self.n_neighbors)
        
        # Step 2: Select landmarks
        if self.landmark_method == 'maxmin':
            landmarks = maxmin_landmarks(X, L)
        elif self.landmark_method == 'random':
            landmarks = random_landmarks(X, L)
        else:
            landmarks = kmeansplusplus_landmarks(X, L)
        
        self.landmarks_ = landmarks
        
        # Step 3: Compute geodesic distances from landmarks
        D_all = np.zeros((L, N))
        for i, landmark in enumerate(landmarks):
            distances = dijkstra(adjacency, indices=landmark)
            D_all[i] = distances
        
        # Handle infinite distances (disconnected components)
        if np.isinf(D_all).any():
            print("Warning: Graph is disconnected. Using finite subgraph.")
            D_all = np.where(np.isinf(D_all), D_all[np.isfinite(D_all)].max() * 2, D_all)
        
        # Step 4: Extract landmark-to-landmark distances
        D_landmarks = D_all[:, landmarks]
        
        # Step 5: Landmark MDS
        self.embedding_, self.eigenvalues_ = landmark_mds(
            D_landmarks, D_all, self.n_components
        )
        
        return self.embedding_

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

import numpy as np
from scipy.sparse.csgraph import dijkstra, connected_components
from scipy.linalg import eigh
 
def robust_landmark_isomap(X, n_neighbors=5, n_landmarks=100, 
                            n_components=2, eigenvalue_tol=1e-10):
    """
    Robust Landmark Isomap with error handling.
    
    Handles:
    - Disconnected graphs
    - Negative eigenvalues
    - Numerical stability issues
    """
    from sklearn.neighbors import kneighbors_graph
    
    N = X.shape[0]
    L = min(n_landmarks, N)
    
    # Build k-NN graph
    adjacency = kneighbors_graph(X, n_neighbors, mode='distance', 
                                  include_self=False)
    adjacency = (adjacency + adjacency.T) / 2  # Symmetrize
    
    # Check connectivity
    n_components_graph, labels = connected_components(
        adjacency, directed=False, return_labels=True
    )
    
    if n_components_graph > 1:
        # Use largest component only
        sizes = np.bincount(labels)
        largest = np.argmax(sizes)
        mask = labels == largest
        
        print(f"Graph has {n_components_graph} components. "
              f"Using largest ({sizes[largest]}/{N} points).")
        
        X_sub = X[mask]
        adjacency_sub = adjacency[mask][:, mask]
        N_sub = X_sub.shape[0]
        
        # Recompute landmarks on connected subset
        landmarks = maxmin_landmarks(X_sub, min(L, N_sub))
    else:
        X_sub = X
        adjacency_sub = adjacency
        N_sub = N
        landmarks = maxmin_landmarks(X, L)
        mask = np.ones(N, dtype=bool)
    
    L = len(landmarks)
    
    # Compute geodesic distances from landmarks
    D_all = np.zeros((L, N_sub))
    for i, landmark in enumerate(landmarks):
        D_all[i] = dijkstra(adjacency_sub, indices=landmark)
    
    # Check for remaining infinite distances
    inf_mask = np.isinf(D_all)
    if inf_mask.any():
        max_finite = D_all[~inf_mask].max() if (~inf_mask).any() else 1.0
        D_all = np.where(inf_mask, max_finite * 2, D_all)
        print(f"Warning: {inf_mask.sum()} infinite distances replaced.")
    
    # Extract landmark-to-landmark distances
    D_landmarks = D_all[:, landmarks]
    
    # Landmark MDS with robust eigenvalue handling
    D_sq = D_landmarks ** 2
    row_mean = D_sq.mean(axis=1, keepdims=True)
    col_mean = D_sq.mean(axis=0, keepdims=True)
    total_mean = D_sq.mean()
    B = -0.5 * (D_sq - row_mean - col_mean + total_mean)
    
    eigenvalues, eigenvectors = eigh(B)
    
    # Sort descending and select top components
    idx = np.argsort(eigenvalues)[::-1][:n_components]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    
    # Handle negative eigenvalues
    n_negative = (eigenvalues < 0).sum()
    if n_negative > 0:
        print(f"Clipping {n_negative} negative eigenvalues to zero.")
        eigenvalues = np.maximum(eigenvalues, 0)
    
    # Nyström extension with numerical stability
    Lambda_inv_sqrt = np.where(
        eigenvalues > eigenvalue_tol,
        1.0 / np.sqrt(eigenvalues),
        0
    )
    
    d_bar = D_sq.mean(axis=1)
    D_all_sq = D_all ** 2
    
    embedding = -0.5 * (D_all_sq.T - d_bar) @ (eigenvectors * Lambda_inv_sqrt)
    
    # Compute residual variance for diagnostics
    residual_variance = compute_residual_variance(D_all, embedding, landmarks)
    
    return embedding, eigenvalues, residual_variance, mask
 
def compute_residual_variance(D_true, embedding, landmark_indices):
    """
    Compute residual variance as diagnostic.
    
    1 - R² between geodesic and embedded distances.
    """
    L, N = D_true.shape
    
    # Compute embedded distances to landmarks
    landmark_embedding = embedding[landmark_indices]
    D_embedded = np.zeros((L, N))
    for i in range(L):
        D_embedded[i] = np.linalg.norm(embedding - landmark_embedding[i], axis=1)
    
    # Compute R²
    D_true_flat = D_true.flatten()
    D_embedded_flat = D_embedded.flatten()
    
    ss_res = np.sum((D_true_flat - D_embedded_flat) ** 2)
    ss_tot = np.sum((D_true_flat - D_true_flat.mean()) ** 2)
    
    r_squared = 1 - ss_res / ss_tot if ss_tot > 0 else 0
    residual_variance = 1 - r_squared
    
    return residual_variance