Coverage for src/jsharpe/sharpe/clustering.py: 100%
77 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 10:09 +0000
1"""Clustering and effective-dimensionality utilities.
3This module groups the routines used to assess the structure of a
4correlation matrix: the effective rank and the optimal number of
5clusters (with silhouette-based quality scoring).
6"""
7# ruff: noqa: N803, N806, S101, TRY003
9import math
10import warnings
12import numpy as np
13import scipy
16def effective_rank(C: np.ndarray) -> float:
17 """Compute the effective rank of a positive semi-definite matrix.
19 The effective rank measures the "effective dimensionality" of a matrix
20 by computing the exponential of the entropy of its normalized eigenvalues.
21 This provides a continuous measure between 1 (perfectly correlated) and
22 n (perfectly uncorrelated/identity matrix).
24 Algorithm:
25 1. Compute eigenvalues (non-negative for PSD matrices)
26 2. Discard zero eigenvalues
27 3. Normalize to form a probability distribution
28 4. Compute entropy H = -Σ p_i log(p_i)
29 5. Return exp(H)
31 Args:
32 C: Positive semi-definite matrix (e.g., correlation matrix).
33 Shape (n, n).
35 Returns:
36 Effective rank, a value in [1, n] where n is the matrix dimension.
38 References:
39 Roy, O. and Vetterli, M. (2007). "The effective rank: a measure of
40 effective dimensionality." EURASIP Journal on Advances in Signal
41 Processing. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.177.2721
43 Example:
44 >>> import numpy as np
45 >>> # Identity matrix has effective rank equal to its dimension
46 >>> abs(effective_rank(np.eye(3)) - 3.0) < 1e-10
47 True
48 >>> # Perfectly correlated matrix has effective rank 1
49 >>> C = np.ones((3, 3))
50 >>> abs(effective_rank(C) - 1.0) < 1e-10
51 True
52 """
53 p = np.linalg.eigvalsh(C)
54 p = p[p > 0]
55 p = p / sum(p)
56 H = np.sum(-p * np.log(p))
57 return math.exp(H)
60def _silhouette_samples(X: np.ndarray, labels: np.ndarray) -> np.ndarray:
61 """Compute silhouette coefficients for each sample.
63 For each sample i:
64 a(i) = mean distance from i to all other samples in the same cluster
65 b(i) = min over other clusters c of mean distance from i to samples in c
66 s(i) = (b(i) - a(i)) / max(a(i), b(i))
68 Args:
69 X: Feature matrix of shape (n_samples, n_features).
70 labels: Cluster labels of shape (n_samples,).
72 Returns:
73 Silhouette coefficients of shape (n_samples,).
74 """
75 n = X.shape[0]
76 unique_labels = np.unique(labels)
77 n_clusters = len(unique_labels)
79 # Pairwise Euclidean distances: dist[i,j] = ||X[i] - X[j]||
80 diff = X[:, np.newaxis, :] - X[np.newaxis, :, :]
81 dist = np.sqrt(np.einsum("ijk,ijk->ij", diff, diff))
83 # For each cluster c, compute mean distance from every sample to cluster c.
84 # For samples IN cluster c, the self-distance (zero) is excluded.
85 label_to_idx = {c: idx for idx, c in enumerate(unique_labels)}
86 cluster_mean_dist = np.zeros((n, n_clusters))
87 for c_idx, c in enumerate(unique_labels):
88 mask = labels == c
89 count = int(mask.sum())
90 dist_to_cluster = dist[:, mask].sum(axis=1)
91 cluster_mean_dist[:, c_idx] = dist_to_cluster / max(count, 1)
92 # For members of cluster c, exclude self-distance (which is 0)
93 in_cluster = np.where(mask)[0]
94 if count > 1:
95 cluster_mean_dist[in_cluster, c_idx] = dist_to_cluster[in_cluster] / (count - 1)
96 else:
97 cluster_mean_dist[in_cluster, c_idx] = 0.0
99 # a(i): mean distance to same cluster
100 own_cluster_idx = np.array([label_to_idx[lbl] for lbl in labels])
101 a = cluster_mean_dist[np.arange(n), own_cluster_idx]
103 # b(i): min mean distance to any other cluster
104 b_mat = cluster_mean_dist.copy()
105 b_mat[np.arange(n), own_cluster_idx] = np.inf
106 b = b_mat.min(axis=1)
107 b = np.where(b == np.inf, 0.0, b)
109 max_ab = np.maximum(a, b)
110 return np.where(max_ab > 0, (b - a) / max_ab, 0.0)
113def _is_correlation_matrix(C: np.ndarray) -> bool:
114 """Return whether ``C`` is a finite correlation matrix.
116 A correlation matrix is a numpy array with all entries in ``[-1, 1]``, ones
117 on the diagonal, and no non-finite values.
119 Args:
120 C: Candidate correlation matrix.
122 Returns:
123 ``True`` if ``C`` satisfies every correlation-matrix property.
124 """
125 return bool(
126 isinstance(C, np.ndarray)
127 and np.all(C >= -1)
128 and np.all(C <= 1)
129 and np.all(np.diag(C) == 1)
130 and np.all(np.isfinite(C))
131 )
134def _validate_correlation_matrix(C: np.ndarray) -> None:
135 """Validate that ``C`` is a finite correlation matrix.
137 Args:
138 C: Candidate correlation matrix.
140 Raises:
141 ValueError: If ``C`` is not a finite correlation matrix (entries in
142 ``[-1, 1]``, ones on the diagonal, all finite).
143 """
144 if not _is_correlation_matrix(C):
145 raise ValueError(
146 "C must be a finite correlation matrix: a square array with all entries in [-1, 1], "
147 "ones on the diagonal, and no non-finite values."
148 )
151def _best_clustering_for_k(D: np.ndarray, k: int, *, retries: int) -> tuple[float, np.ndarray | None]:
152 """Run k-means ``retries`` times for a fixed ``k`` and return the best solution.
154 Degenerate runs (fewer than ``k`` non-empty clusters, or silhouette scores
155 with zero spread) are skipped. Quality is the mean silhouette score divided
156 by its standard deviation.
158 Args:
159 D: Distance matrix used as the feature matrix for k-means.
160 k: Number of clusters to fit.
161 retries: Number of k-means restarts to reduce sensitivity to the random
162 initialisation.
164 Returns:
165 Tuple of (quality, labels) for the best valid run, or ``(-inf, None)``
166 if no run produced a valid clustering.
167 """
168 best_quality = -np.inf
169 best_labels: np.ndarray | None = None
170 for _ in range(retries):
171 with warnings.catch_warnings():
172 warnings.simplefilter("ignore")
173 _centroids, labels = scipy.cluster.vq.kmeans2(D, k, minit="points", iter=300)
174 # Skip degenerate solutions with empty clusters
175 if len(np.unique(labels)) < k:
176 continue
177 silhouette_vals = _silhouette_samples(D, labels)
178 std = silhouette_vals.std()
179 if std == 0:
180 continue
181 q = float(silhouette_vals.mean() / std)
182 if q > best_quality:
183 best_quality = q
184 best_labels = labels.copy()
185 return best_quality, best_labels
188def _select_best_k(qualities: dict[int, float], best_labels: dict[int, np.ndarray]) -> int:
189 """Return the ``k`` with the highest quality among those with a valid clustering.
191 Args:
192 qualities: Mapping of ``k`` to its quality score (``-inf`` for degenerate ``k``).
193 best_labels: Mapping of ``k`` to its best label assignment; keys mark the
194 values of ``k`` for which a valid clustering was found.
196 Returns:
197 The number of clusters ``k`` maximising the quality score.
199 Raises:
200 RuntimeError: If no ``k`` produced a valid clustering.
201 """
202 valid_k = {k: q for k, q in qualities.items() if k in best_labels}
203 if not valid_k:
204 raise RuntimeError("No valid clustering solution found; try increasing retries or reducing max_clusters.")
205 return max(valid_k, key=lambda x: valid_k[x])
208def number_of_clusters(
209 C: np.ndarray,
210 *,
211 retries: int = 10,
212 max_clusters: int = 100,
213) -> tuple[int, dict[int, float], np.ndarray]:
214 """Compute the optimal number of clusters from a correlation matrix.
216 Implements the algorithm from section 8.1 of Lopez de Prado (2018):
217 1. Convert the correlation matrix into a distance matrix.
218 2. Using the columns of the distance matrix as features, run the
219 k-means algorithm for each k and compute the quality of the
220 clustering.
221 3. Return the clustering with the highest quality.
223 Quality is defined as the mean of the silhouette scores divided by
224 their standard deviation.
226 Args:
227 C: Correlation matrix. Must be square, symmetric, finite, with ones
228 on the diagonal and all entries in [-1, 1].
229 retries: Number of times to run k-means for each k to reduce the
230 impact of random initialisation. Default 10.
231 max_clusters: Maximum number of clusters to evaluate. Capped at
232 ``C.shape[0] - 1``. Default 100.
234 Returns:
235 Tuple of (n_clusters, qualities, labels):
236 - n_clusters: Optimal number of clusters.
237 - qualities: Dict mapping k to its quality score.
238 - labels: Cluster assignment for each observation (shape (n,)).
240 Raises:
241 ValueError: If ``C`` is not a finite correlation matrix.
242 RuntimeError: If no ``k`` produced a valid clustering.
244 References:
245 Lopez de Prado, M. (2018). "Detection of false investment strategies
246 using unsupervised learning methods." SSRN 3167017.
247 https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3167017
249 Example:
250 >>> from jsharpe.sharpe.generators import get_random_correlation_matrix
251 >>> np.random.seed(42)
252 >>> C, _, _ = get_random_correlation_matrix(
253 ... number_of_trials=20, effective_number_of_trials=4
254 ... )
255 >>> n, qualities, labels = number_of_clusters(C, retries=3, max_clusters=8)
256 >>> 2 <= n <= 8
257 True
258 >>> labels.shape
259 (20,)
260 """
261 _validate_correlation_matrix(C)
263 max_clusters = min(max_clusters, C.shape[0] - 1)
265 # Convert correlations to distances. C was validated above to be a finite
266 # correlation matrix (entries in [-1, 1]), so (1 - C) / 2 lies in [0, 1] and
267 # its square root is always finite: this assert is a pure internal invariant
268 # a caller cannot trigger, not a public-API guard.
269 D = np.sqrt((1 - C) / 2)
270 assert np.all(np.isfinite(D))
272 qualities: dict[int, float] = {}
273 best_labels: dict[int, np.ndarray] = {}
274 for k in range(2, max_clusters + 1):
275 quality, labels = _best_clustering_for_k(D, k, retries=retries)
276 qualities[k] = quality
277 if labels is not None:
278 best_labels[k] = labels
280 # Select the best k among those for which a valid solution was found
281 best_k = _select_best_k(qualities, best_labels)
282 return best_k, qualities, best_labels[best_k]