Spectral Clustering

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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313

import numpy as np
from scipy.linalg import eigh
from scipy.sparse.linalg import eigsh
from scipy.sparse import diags, csr_matrix
from scipy.spatial.distance import cdist
from sklearn.neighbors import NearestNeighbors
from sklearn.cluster import KMeans
from sklearn.preprocessing import normalize
from typing import Optional, Literal, Tuple
import warnings
 
class SpectralClusteringComplete:
    """
    Complete implementation of spectral clustering with all options.
    
    This implementation makes every choice explicit and provides
    diagnostics for understanding clustering quality.
    """
    
    def __init__(
        self,
        n_clusters: int,
        affinity: Literal['rbf', 'knn', 'precomputed'] = 'rbf',
        gamma: Optional[float] = None,
        n_neighbors: int = 10,
        laplacian: Literal['unnormalized', 'symmetric', 'random_walk'] = 'symmetric',
        row_normalize: bool = True,
        n_init: int = 10,
        random_state: Optional[int] = None
    ):
        """
        Initialize spectral clustering.
        
        Parameters:
        -----------
        n_clusters : int
            Number of clusters k
        affinity : str
            How to construct similarity graph:
            - 'rbf': Gaussian kernel
            - 'knn': K-nearest neighbors
            - 'precomputed': W is provided directly
        gamma : float, optional
            RBF kernel parameter (= 1/(2σ²)). If None, estimated from data.
        n_neighbors : int
            Number of neighbors for KNN affinity
        laplacian : str
            Which Laplacian to use:
            - 'unnormalized': L = D - W (RatioCut)
            - 'symmetric': L_sym = I - D^{-1/2}WD^{-1/2} (NCut)
            - 'random_walk': L_rw = I - D^{-1}W (NCut)
        row_normalize : bool
            Whether to row-normalize eigenvector embedding (NJW style)
        n_init : int
            Number of K-Means initializations
        random_state : int, optional
            Random seed for reproducibility
        """
        self.n_clusters = n_clusters
        self.affinity = affinity
        self.gamma = gamma
        self.n_neighbors = n_neighbors
        self.laplacian = laplacian
        self.row_normalize = row_normalize
        self.n_init = n_init
        self.random_state = random_state
        
        # Stored for diagnostics
        self.W_ = None
        self.eigenvalues_ = None
        self.eigenvectors_ = None
        self.embedding_ = None
        
    def _construct_similarity_graph(self, X: np.ndarray) -> np.ndarray:
        """Construct similarity matrix W from data X."""
        n = X.shape[0]
        
        if self.affinity == 'precomputed':
            # X is already the similarity matrix
            W = X.copy()
            np.fill_diagonal(W, 0)
            return W
        
        elif self.affinity == 'rbf':
            # Gaussian kernel
            distances = cdist(X, X, metric='euclidean')
            
            if self.gamma is None:
                # Estimate gamma from data
                # Common heuristic: median of pairwise distances
                nonzero_dists = distances[distances > 0]
                sigma = np.median(nonzero_dists)
                gamma = 1.0 / (2 * sigma**2)
            else:
                gamma = self.gamma
            
            W = np.exp(-gamma * distances**2)
            np.fill_diagonal(W, 0)
            return W
        
        elif self.affinity == 'knn':
            # K-nearest neighbors graph
            nn = NearestNeighbors(n_neighbors=self.n_neighbors + 1)
            nn.fit(X)
            
            A = nn.kneighbors_graph(mode='connectivity')
            A = A.toarray()
            np.fill_diagonal(A, 0)
            
            # Symmetrize
            W = 0.5 * (A + A.T)
            return W
        
        else:
            raise ValueError(f"Unknown affinity: {self.affinity}")
    
    def _compute_laplacian(self, W: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        Compute the specified Laplacian matrix.
        
        Returns:
        --------
        L : ndarray
            Laplacian matrix
        d : ndarray
            Degree vector
        """
        n = W.shape[0]
        d = np.sum(W, axis=1)
        
        # Handle isolated nodes (zero degree)
        isolated = d == 0
        if np.any(isolated):
            warnings.warn(f"{np.sum(isolated)} isolated nodes detected. "
                         "Consider adjusting affinity parameters.")
            d[isolated] = 1.0  # Prevent division by zero
        
        if self.laplacian == 'unnormalized':
            L = np.diag(d) - W
            
        elif self.laplacian == 'symmetric':
            d_inv_sqrt = 1.0 / np.sqrt(d)
            D_inv_sqrt = np.diag(d_inv_sqrt)
            L = np.eye(n) - D_inv_sqrt @ W @ D_inv_sqrt
            
        elif self.laplacian == 'random_walk':
            d_inv = 1.0 / d
            D_inv = np.diag(d_inv)
            L = np.eye(n) - D_inv @ W
            
        else:
            raise ValueError(f"Unknown Laplacian type: {self.laplacian}")
        
        return L, d
    
    def _eigendecomposition(self, L: np.ndarray, k: int) -> Tuple[np.ndarray, np.ndarray]:
        """
        Compute k smallest eigenvalues/eigenvectors.
        
        For small matrices, uses dense eigendecomposition.
        For large sparse matrices, use scipy.sparse.linalg.eigsh.
        """
        n = L.shape[0]
        
        if n <= 1000:
            # Dense eigendecomposition
            eigenvalues, eigenvectors = eigh(L)
            # Take first k
            eigenvalues = eigenvalues[:k]
            eigenvectors = eigenvectors[:, :k]
        else:
            # Sparse eigendecomposition (for larger matrices)
            # eigsh finds smallest algebraic eigenvalues
            L_sparse = csr_matrix(L)
            eigenvalues, eigenvectors = eigsh(L_sparse, k=k, which='SM')
            
            # Sort by eigenvalue
            idx = np.argsort(eigenvalues)
            eigenvalues = eigenvalues[idx]
            eigenvectors = eigenvectors[:, idx]
        
        return eigenvalues, eigenvectors
    
    def fit_predict(self, X: np.ndarray) -> np.ndarray:
        """
        Perform spectral clustering.
        
        Parameters:
        -----------
        X : ndarray, shape (n_samples, n_features)
            Data to cluster (or precomputed affinity matrix if affinity='precomputed')
        
        Returns:
        --------
        labels : ndarray, shape (n_samples,)
            Cluster assignments
        """
        # Step 1: Construct similarity graph
        self.W_ = self._construct_similarity_graph(X)
        
        # Step 2: Compute Laplacian
        L, self.degrees_ = self._compute_laplacian(self.W_)
        
        # Step 3: Eigendecomposition
        self.eigenvalues_, self.eigenvectors_ = self._eigendecomposition(
            L, self.n_clusters
        )
        
        # Step 4: Form embedding
        U = self.eigenvectors_.copy()
        
        # Step 5: Row normalization (optional but recommended)
        if self.row_normalize:
            row_norms = np.linalg.norm(U, axis=1, keepdims=True)
            # Handle zero rows
            zero_rows = row_norms.flatten() < 1e-10
            row_norms[zero_rows] = 1.0
            U = U / row_norms
        
        self.embedding_ = U
        
        # Step 6: K-Means clustering
        kmeans = KMeans(
            n_clusters=self.n_clusters,
            n_init=self.n_init,
            random_state=self.random_state
        )
        labels = kmeans.fit_predict(U)
        
        self.labels_ = labels
        self.kmeans_ = kmeans
        
        return labels
    
    def get_diagnostics(self) -> dict:
        """
        Return diagnostic information for assessing clustering quality.
        """
        if self.eigenvalues_ is None:
            raise ValueError("Must call fit_predict first")
        
        # Compute eigenvalue gaps
        n_eig = len(self.eigenvalues_)
        if n_eig > 1:
            gaps = np.diff(self.eigenvalues_)
        else:
            gaps = np.array([])
        
        # Count near-zero eigenvalues
        near_zero = np.sum(self.eigenvalues_ < 1e-10)
        
        # Estimate number of disconnected components
        num_components = near_zero
        
        return {
            'eigenvalues': self.eigenvalues_,
            'eigenvalue_gaps': gaps,
            'spectral_gap': self.eigenvalues_[-1] if n_eig > 0 else None,
            'num_near_zero_eigenvalues': near_zero,
            'estimated_components': num_components,
            'graph_density': np.sum(self.W_ > 0) / (self.W_.shape[0]**2),
            'degree_stats': {
                'min': self.degrees_.min(),
                'max': self.degrees_.max(),
                'mean': self.degrees_.mean(),
                'std': self.degrees_.std()
            }
        }
 
 
def demonstrate_complete_algorithm():
    """
    Demonstrate the complete spectral clustering algorithm.
    """
    from sklearn.datasets import make_moons, make_circles
    from sklearn.metrics import adjusted_rand_score, silhouette_score
    
    print("Complete Spectral Clustering Algorithm")
    print("=" * 60)
    
    # Test on different datasets
    datasets = [
        ("Two Moons", *make_moons(300, noise=0.05, random_state=42), 2),
        ("Concentric Circles", *make_circles(300, noise=0.05, factor=0.5, random_state=42), 2),
    ]
    
    for name, X, y_true, k in datasets:
        print(f"\n--- {name} ---")
        
        sc = SpectralClusteringComplete(
            n_clusters=k,
            affinity='rbf',
            laplacian='symmetric',
            row_normalize=True,
            random_state=42
        )
        
        labels = sc.fit_predict(X)
        
        ari = adjusted_rand_score(y_true, labels)
        sil = silhouette_score(X, labels)
        
        print(f"ARI: {ari:.4f}, Silhouette: {sil:.4f}")
        
        # Diagnostics
        diag = sc.get_diagnostics()
        print(f"Eigenvalues: {diag['eigenvalues'].round(4)}")
        print(f"Graph density: {diag['graph_density']:.4f}")
        print(f"Degree range: [{diag['degree_stats']['min']:.2f}, {diag['degree_stats']['max']:.2f}]")
 
 
if __name__ == "__main__":
    demonstrate_complete_algorithm()

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
131
132
133
134
135
136
137
138
139
140
141

import numpy as np
from sklearn.cluster import SpectralClustering
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score
from scipy.spatial.distance import cdist
import warnings
 
def analyze_sigma_sensitivity():
    """
    Analyze sensitivity to the Gaussian kernel bandwidth parameter.
    """
    print("Sensitivity to Kernel Bandwidth (σ)")
    print("=" * 60)
    
    X, y_true = make_moons(300, noise=0.05, random_state=42)
    
    # Compute distance statistics for reference
    distances = cdist(X, X)
    nonzero_dists = distances[distances > 0]
    
    print(f"\nDistance statistics:")
    print(f"  Min: {nonzero_dists.min():.4f}")
    print(f"  Median: {np.median(nonzero_dists):.4f}")
    print(f"  Mean: {nonzero_dists.mean():.4f}")
    print(f"  Max: {nonzero_dists.max():.4f}")
    
    print(f"\nGamma = 1/(2σ²), so larger gamma = smaller σ")
    print("\nGamma    |   σ      | ARI    | Interpretation")
    print("-" * 60)
    
    for gamma in [0.01, 0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 50.0, 100.0]:
        sigma = np.sqrt(1 / (2 * gamma))
        
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            sc = SpectralClustering(
                n_clusters=2, affinity='rbf', gamma=gamma,
                random_state=42, n_init=10
            )
            labels = sc.fit_predict(X)
        
        ari = adjusted_rand_score(y_true, labels)
        
        if gamma < 0.1:
            interp = "σ too large, graph nearly complete"
        elif gamma > 20:
            interp = "σ too small, graph sparse/disconnected"
        elif 0.5 <= gamma <= 5:
            interp = "Good range"
        else:
            interp = "Marginal"
        
        print(f"{gamma:8.2f} | {sigma:8.4f} | {ari:.4f} | {interp}")
 
 
def analyze_k_neighbor_sensitivity():
    """
    Analyze sensitivity to the number of neighbors in KNN graph.
    """
    print("\n" + "=" * 60)
    print("Sensitivity to Number of Neighbors (k_nn)")
    print("=" * 60)
    
    X, y_true = make_moons(300, noise=0.05, random_state=42)
    n = len(X)
    
    print(f"\nn = {n}, recommended k_nn range: {int(np.log(n))}-{int(np.sqrt(n))}")
    print("\nk_nn     | ARI    | Notes")
    print("-" * 45)
    
    for k_nn in [2, 3, 5, 10, 15, 20, 30, 50, 100]:
        if k_nn >= n:
            continue
            
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            sc = SpectralClustering(
                n_clusters=2, affinity='nearest_neighbors', 
                n_neighbors=k_nn, random_state=42, n_init=10
            )
            labels = sc.fit_predict(X)
        
        ari = adjusted_rand_score(y_true, labels)
        
        if k_nn < 5:
            notes = "May be too sparse"
        elif k_nn > 50:
            notes = "May lose local structure"
        else:
            notes = "Reasonable range"
        
        print(f"{k_nn:8d} | {ari:.4f} | {notes}")
 
 
def demonstrate_eigengap_for_k():
    """
    Demonstrate using the eigengap to choose k.
    """
    from sklearn.datasets import make_blobs
    from scipy.linalg import eigh
    
    print("\n" + "=" * 60)
    print("Choosing k via Eigengap Heuristic")
    print("=" * 60)
    
    # Generate data with known k
    for true_k in [2, 4, 6]:
        X, _ = make_blobs(n_samples=200, centers=true_k, 
                         cluster_std=0.5, random_state=42)
        
        # Build similarity graph
        distances = cdist(X, X)
        sigma = np.median(distances[distances > 0])
        W = np.exp(-distances**2 / (2 * sigma**2))
        np.fill_diagonal(W, 0)
        
        # Compute normalized Laplacian
        d = np.sum(W, axis=1)
        d_inv_sqrt = 1.0 / np.sqrt(np.maximum(d, 1e-10))
        L_sym = np.eye(len(X)) - np.diag(d_inv_sqrt) @ W @ np.diag(d_inv_sqrt)
        
        eigenvalues, _ = eigh(L_sym)
        eigenvalues = eigenvalues[:12]
        
        # Find gaps
        gaps = np.diff(eigenvalues)
        
        # Largest gap after the first few
        suggested_k = np.argmax(gaps[1:8]) + 2  # Skip first, add 2 for 1-indexing + shift
        
        print(f"\nTrue k = {true_k}:")
        print(f"  Eigenvalues: {eigenvalues[:8].round(4)}")
        print(f"  Gaps: {gaps[:7].round(4)}")
        print(f"  Suggested k (largest gap): {suggested_k}")
        print(f"  Correct: {'✓' if suggested_k == true_k else '✗'}")
 
 
if __name__ == "__main__":
    analyze_sigma_sensitivity()
    analyze_k_neighbor_sensitivity()
    demonstrate_eigengap_for_k()

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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

import numpy as np
from scipy.linalg import eigh
from scipy.spatial.distance import cdist
from sklearn.cluster import SpectralClustering
from sklearn.metrics import silhouette_score
import warnings
 
def diagnose_spectral_clustering(X, n_clusters, affinity='rbf', gamma=None):
    """
    Diagnose potential issues with spectral clustering.
    
    Returns diagnostic information to help identify failure modes.
    """
    n = X.shape[0]
    
    # Build similarity graph
    if affinity == 'rbf':
        distances = cdist(X, X)
        if gamma is None:
            sigma = np.median(distances[distances > 0])
            gamma = 1.0 / (2 * sigma**2)
        else:
            sigma = np.sqrt(1 / (2 * gamma))
        
        W = np.exp(-gamma * distances**2)
        np.fill_diagonal(W, 0)
    else:
        raise ValueError("Only 'rbf' affinity supported in diagnostics")
    
    # Compute Laplacian
    d = np.sum(W, axis=1)
    d_inv_sqrt = np.zeros(n)
    d_inv_sqrt[d > 0] = 1.0 / np.sqrt(d[d > 0])
    L_sym = np.eye(n) - np.diag(d_inv_sqrt) @ W @ np.diag(d_inv_sqrt)
    
    # Eigendecomposition
    eigenvalues, eigenvectors = eigh(L_sym)
    
    # Diagnostics
    diagnostics = {
        'n_samples': n,
        'n_clusters_requested': n_clusters,
        'sigma': sigma,
        'gamma': gamma,
    }
    
    # Connectivity check
    near_zero_threshold = 1e-10
    n_zero_eigenvalues = np.sum(eigenvalues < near_zero_threshold)
    diagnostics['n_connected_components'] = n_zero_eigenvalues
    
    if n_zero_eigenvalues > 1:
        diagnostics['warning_disconnected'] = True
        diagnostics['disconnected_message'] = (
            f"Graph has {n_zero_eigenvalues} connected components. "
            "Consider increasing sigma/gamma or using more neighbors."
        )
    else:
        diagnostics['warning_disconnected'] = False
    
    # Eigengap analysis
    eigenvalues_show = eigenvalues[:min(10, n)]
    gaps = np.diff(eigenvalues_show)
    diagnostics['first_10_eigenvalues'] = eigenvalues_show
    diagnostics['eigenvalue_gaps'] = gaps
    
    # Check for clear gap
    if n_clusters < len(gaps):
        gap_at_k = gaps[n_clusters - 1] if n_clusters > 1 else gaps[0]
        median_gap = np.median(gaps[1:])  # Exclude first gap
        
        if gap_at_k < median_gap:
            diagnostics['warning_weak_gap'] = True
            diagnostics['weak_gap_message'] = (
                f"Gap at k={n_clusters} ({gap_at_k:.4f}) is below median gap ({median_gap:.4f}). "
                "Cluster structure may not be clear."
            )
        else:
            diagnostics['warning_weak_gap'] = False
    
    # Degree variance check
    degree_cv = np.std(d) / np.mean(d)  # Coefficient of variation
    diagnostics['degree_cv'] = degree_cv
    
    if degree_cv > 2.0:
        diagnostics['warning_high_degree_variance'] = True
        diagnostics['degree_variance_message'] = (
            f"High degree variance (CV={degree_cv:.2f}). "
            "Ensure using normalized Laplacian."
        )
    else:
        diagnostics['warning_high_degree_variance'] = False
    
    # Run clustering and get silhouette
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        sc = SpectralClustering(
            n_clusters=n_clusters, affinity='rbf', gamma=gamma,
            random_state=42, n_init=10
        )
        labels = sc.fit_predict(X)
    
    sil = silhouette_score(X, labels)
    diagnostics['silhouette_score'] = sil
    
    if sil < 0.2:
        diagnostics['warning_low_silhouette'] = True
        diagnostics['silhouette_message'] = (
            f"Low silhouette score ({sil:.3f}). "
            "Clusters may be overlapping or parameters may be wrong."
        )
    else:
        diagnostics['warning_low_silhouette'] = False
    
    diagnostics['labels'] = labels
    
    return diagnostics
 
 
def print_diagnostics(diag):
    """Pretty-print diagnostic results."""
    print("\nSpectral Clustering Diagnostics")
    print("=" * 55)
    print(f"Samples: {diag['n_samples']}, Clusters: {diag['n_clusters_requested']}")
    print(f"Sigma: {diag['sigma']:.4f}, Gamma: {diag['gamma']:.4f}")
    print(f"\nEigenvalues: {diag['first_10_eigenvalues'][:6].round(4)}")
    print(f"Gaps: {diag['eigenvalue_gaps'][:5].round(4)}")
    print(f"Connected components: {diag['n_connected_components']}")
    print(f"Degree CV: {diag['degree_cv']:.3f}")
    print(f"Silhouette: {diag['silhouette_score']:.3f}")
    
    print("\nWarnings:")
    warnings_found = False
    for key in ['disconnected', 'weak_gap', 'high_degree_variance', 'low_silhouette']:
        if diag.get(f'warning_{key}', False):
            print(f"  ⚠ {diag[f'{key}_message']}")
            warnings_found = True
    
    if not warnings_found:
        print("  ✓ No issues detected")
 
 
def demonstrate_diagnostics():
    """Demonstrate diagnostics on various scenarios."""
    from sklearn.datasets import make_moons, make_blobs
    
    print("Failure Mode Diagnostics Demonstration")
    print("=" * 60)
    
    # Good case: Two moons
    print("\n--- Case 1: Two Moons (should work well) ---")
    X, _ = make_moons(200, noise=0.05, random_state=42)
    diag = diagnose_spectral_clustering(X, n_clusters=2)
    print_diagnostics(diag)
    
    # Bad case: Wrong k
    print("\n--- Case 2: Two Moons with k=4 (wrong k) ---")
    diag = diagnose_spectral_clustering(X, n_clusters=4)
    print_diagnostics(diag)
    
    # Bad case: Overlapping clusters
    print("\n--- Case 3: Overlapping Blobs (hard to separate) ---")
    X_overlap, _ = make_blobs(200, centers=2, cluster_std=2.0, random_state=42)
    diag = diagnose_spectral_clustering(X_overlap, n_clusters=2)
    print_diagnostics(diag)
 
 
if __name__ == "__main__":
    demonstrate_diagnostics()