Improving Representation Learning for Histopathologic Images with Cluster Constraints

2.1. Self-supervised learning

Self-supervised learning (SSL) methods have recently demonstrated effectiveness for computer vision tasks by learning representations without reliance on manual labels. SSL techniques leverage the intrinsic structure and consistency of the data itself as a supervisory signal. Several paradigms have arisen, including contrastive-based, knowledge-distillation-based, clustering-based, and information maximization-based approaches. Typically, these techniques function by generating augmented pairs of views from a single data instance and directing the model to produce similar outputs for each view.

Contrastive learning stands out as a key self-supervised approach in representation learning, with the goal of deriving informative and concise representations from unstructured data. A slew of methods grounded in contrastive learning have been proposed, including Contrastive Predictive Coding (CPC) [35], SimCLR [11], MoCo [19], and NNCLR [17]. CPC, recognized for its widespread application, employs an autoregressive model to predict future observations based on past observations, rendering it particularly apt for sequential data. MoCo and SimCLR, two other popular contrastive learning techniques for instance discrimination, generate positive pairs by utilizing two different views (augmentations) of the same image and negative pairs by pairing augmentations of different data points. The principle behind contrastive learning is to distance negative pairs while converging positive pairs. However, achieving optimal performance with these methods often necessitates a plethora of negative pairs. MoCo addresses this issue by implementing momentum encoders and a memory bank mechanism, while SimCLR capitalizes on large batch sizes for negative pair comparisons. In NNCLR, the model learns representations by minimizing the distance between an anchor and its nearest neighbor in the momentum encoder’s output space while maximizing the distance to other negative samples. This method streamlines the utilization of negative samples in a batch, curtailing the need for large batch sizes and memory banks while maintaining competitive performance compared to other contrastive learning methods.

Knowledge-distillation-based methodologies, such as BYOL [18] and SimSiam [12], aim to enhance performance with smaller batch sizes and without the need for negative samples. In stark contrast to their contrastive counterparts, these non-contrastive techniques employ only positive pairs, eliminating the demand for large batch sizes or a memory bank mechanism. BYOL stands out with its momentum update mechanism, which renders negative pairs unnecessary. It establishes a target network by applying an exponential moving average to the online network’s weights. This “moving target” offers the online network a stable benchmark to aim for during training, pushing the network away from trivial solutions and encouraging richer representations. SimSiam, influenced by BYOL, streamlines the process by forgoing the moving target. Instead, it adopts a symmetric architecture, where dual networks reciprocally predict each other’s outputs. To ensure the representations are non-trivial, SimSiam utilizes a stop-gradient operation. Nonetheless, despite their batch efficiency, both BYOL and SimSiam remain susceptible to collapsing into trivial solutions.

Clustering-based SSL methods, such as DeepCluster [4] and Prototypical Contrastive Learning (PCL) [31], offer an innovative angle to representation learning by capitalizing on iterative pseudo-labeling to cluster the learned representations. DeepCluster updates its network parameters according to pseudo-labels produced by the K-Means clustering algorithm. This aligns the network more closely with the inherent data distribution. Expanding on the foundation laid by MoCo, PCL integrates the ProtoNCE loss and K-Means clustering, aiming to refine image embeddings by nudging them closer to their respective prototypes by optimizing the ProtoNCE loss function. SwAV [5], another clustering-based SSL method, shares similarities with SimSiam architecture but differentiates itself with a swap prediction mechanism. Specifically, SwAV aligns the cluster assignments of one augmentation with the representations of another version of the same image. These assignments are fine-tuned using the Sinkhorn algorithm, ensuring a balanced distribution. In embracing this strategy, SwAV learns invariant features that encapsulate crucial semantic information.

Information maximization-based SSL approaches strive to learn representations by maximizing the invariance of corresponding features while minimizing the covariance between different features. Noteworthy methods in this realm include VICReg [2] and Barlow Twins [45]. VICReg specifically enforces feature representation invariance by amplifying their variance and curtailing the covariance between different features. This strategy ensures that the learned representations are not only informative but also capture the core attributes of the data. Contrasting with VICReg, which explicitly maximizes the variance of individual features, Barlow Twins centers its focus on minimizing cross-correlations across feature dimensions. It achieves this by mitigating the cross-correlation between the outputs of twin networks, each processing a distinct augmentation of the same image while concurrently accentuating the invariance of matching features.

2.2. SSL in Pathology

The advent of whole-slide scanners has enabled the digitization of histopathological slides, gradually transforming the field of anatomical pathology into a data-abundant domain. Recognizing this, SSL techniques are being increasingly employed in computational pathology to take advantage of the abundance of unlabeled data. In conjunction with smaller labeled datasets, these techniques promise to elevate diagnostic precision and bolster predictive modeling.

In the realm of computational pathology, SSL is gradually gaining adoption as a means to tackle challenges such as acquiring annotations for pathology slides, managing high-resolution images, and addressing the substantial variability in their features. Numerous studies have leveraged SSL for extracting features from WSIs and have utilized these features to achieve promising results in downstream histopathology image analysis tasks [10, 16, 26, 30]. With the escalating adoption of SSL in computational pathology, there is a growing necessity to determine the applicability of general self-supervised methods to histopathology images. A recent benchmarking study [27] gauged various SSL methods across diverse pathology datasets for various downstream tasks, such as classification and nuclei instance segmentation tasks. Their results indicate that SSL can considerably uplift the performance on downstream tasks on histopathology images compared to ImageNet pre-trained and supervised models, especially when labeled data is scarce.

3.1. Preliminaries: SimSiam

Self-supervised visual representation learning is a method for learning an embedding function that maps an input image $x$ to a representation. This is typically achieved by using a similarity measure designed to enforce similarity between augmented views. Starting with sets of data transformations, ${𝒯}_1$ and ${𝒯}_2$, we randomly sample transformations $t_1,t_2~{𝒯}_1,{𝒯}_2$ and produce augmented views $x_1=t_1\left(x\right)$ and $x_2=t_2\left(x\right)$. An encoder $f$ is used to produce representations $y_1=f\left(x_1\right)$ and $y_2=f\left(x_2\right)$, which are then fed to a projector $h$ to produce projections $z_1=h\left(y_1\right)$ and $z_2=h\left(y_2\right)$. As in SimSiam, we pass $z_1$ through a predictor $g$ to produce the prediction $p_1=g\left(z_1\right)$. Additionally, we swap the views and produce a symmetric loss as follows:

where

$\parallel \cdot \parallel$ is the ${\ell }_2-$ norm, and $\text{sg}(\cdot )$ is the stop-gradient operation to prevent collapse.

3.2. CluSiam: Cluster-Constrained SSL

We build upon the SimSiam architecture by adding a cluster assigner $q$ that operates on the projections produced by $h$. We use the outputs of $h$ as the inputs to $q$ because the batch normalization layers $\left(ℬ𝒩\right)$ in $h$ stabilize the distribution of $q$’s inputs. We also introduce $ℬ𝒩$ in $q$ because controlling the input scale is crucial for generating cluster assignment probabilities using softmax. Given batches of views $X_1,X_2$, we produce projections $Z_1=\left\{h\left(f\left(x_i^{\left(1\right)}\right)\right)\right. : \left.x_i^{\left(1\right)}\in X_1\right\},Z_2=\left\{h\left(f\left(x_i^{\left(2\right)}\right)\right):x_i^{\left(2\right)}\in X_2\right\}$. We then concatenate these projections to obtain $Z=\text{concat}\left(Z_1,Z_2\right)$. The cluster assigner $q$ maps the concatenated projections $Z$ to cluster representations $A$ defined as:

where $A_{ij}$ is the element at the $i$-th row and $j$-th column of $A,\tau$ is the temperature, and $K$ is the exploration space that represents the maximum number of clusters allowed during the training. This operation is applied along each row, meaning that for every data point $i$, the sum of $A_{ij}$ over all clusters, $j$’s, equals 1. Finally, we map the cluster representation to the clusters as follows:

Notably, $\text{argmax}\left(A\right)$ is not differentiable, so there is no real gradient for this operation. We approximate the gradient similar to the straight-through estimator and just copy gradients from $A$ to $\text{argmax}\left(A\right),{\nabla }_{\text{argmax}\left(A\right)}C={\nabla }_{A}C$, making it possible for backpropagation. $\text{sg}(\cdot )$ denotes the stop-gradient operation, $\text{nonzero}(·)$ filters out vectors along the row dimension that are all zeros, and $k$ represents the count of non-zero centroids. The dimension $D$ corresponds to the feature dimension and is consistent with the dimensions of $p$ and $z$.

In our clustering module, unlike common clustering methods that use inter-class similarity or other SSL techniques that focus on the invariance between two different views, we do not impose any restrictions on inter-class similarity or the assignments between two different views, $a_i$ and $a_{n+i}$. Our cluster loss is solely defined by inter-cluster separation. This separation is defined as:

Importantly, all the tensors in (5) are ${\ell }_2$ normalized, so equation (5) can be interpreted as representing the average cosine similarities between clusters. Furthermore, the stop-gradient operation is applied to all elements in (4) to prevent the cluster assigner from collapsing to a trivial solution where all samples are assigned to the same set or cluster. We can view $A=a_1,\cdots ,a_{2n}$ as a latent variable, with our goal being to minimize $ℒ(A,Z)$. This optimization problemcan be solved by an alternating algorithm that fixes one set of variables and solves for the other set.

The intra-class similarity term was not introduced in this design because we neither use a contrastive formulation loss function like SimCLR nor an additional projection head to introduce the knowledge distillation architecture like BYOL. Directly optimizing for high intra-class similarity is prone to collapsing all samples into the same representation, which is a trivial solution and does not capture any meaningful information.

However, our cluster learning is still prone to collapse due to the presence of the softmax and argmax functions. The hard cluster assignment of the softmax and argmax functions limits the model to updating only the maximum scoring cluster during training. This does not encourage exploration of different combinations of cluster assignments, and other potential cluster assignments are left out of the backpropagation updates. As a result, the model is prone to getting stuck in a trivial single-cluster solution during training, a phenomenon known as “collapsing”. This is especially likely as the number of valid clusters decreases. Using softmax and argmax based cluster assignment can exacerbate this problem, leading to the rapid collapse of a single cluster. This prevents the model from learning the intended cluster structures and, in turn, leads to meaningless representations since the cluster loss cannot be effectively optimized with only a single active cluster.

The task of cluster assignment can be viewed as a decision-making process in which the cluster assigner determines which samples belong in the same cluster. Striking a balance between exploration and exploitation is crucial in this decision-making process. Exploration refers to trialing different actions to learn more about the environment and their associated losses, while exploitation refers to choosing the action currently known to possess the lowest expected loss. To prevent the collapse of clustering due to continual updates only to the neuron with the highest probability value, we introduce randomness into the decision-making process by adding Gumbel noise [23]. Gumbel noise is a random variable sampled from a Gumbel distribution. It proves useful in discrete action spaces, where a model must choose between a finite set of actions, as in our cluster assignment task. By replacing (3) with (6), the cluster assigner can explore different cluster combinations based on their probabilities, thereby learning more about different cluster combinations and their associated losses.

where $G$ represents the Gumbel distribution.

Therefore, CluSiam can be trained effectively using the composite loss function, as shown in equation (7). In this equation, $\alpha$ serves as a hyperparameter, adjusting the magnitude of the weights. In our implementation, we set $\alpha =0.5$ by default.

4.1. Clustering

In the clustering task, we evaluated the performance of clustering algorithms using the Rand Index (RI). The performance was compared using different representations and clustering algorithms (Table 1). Importantly, our cluster assigner’s output is different from traditional methods like K-Means. Unlike K-Means, which provides a hard cluster assignment for each patch, our assigner outputs a K-dimensional vector. This structure allows for more flexibility in generating cluster assignments beyond the typical use of argmax. Such probabilistic outputs grant our method greater adaptability in clustering, especially when contrasted with the rigid assignments derived from K-Means. We visualized the cluster assignments of two WSIs generated simply using argmax, alongside their respective ground truths, in Figure 3.

4.2. Classification

For the classification task, we evaluated model performance using accuracy and Area Under the ROC Curve (AUC) metrics. We aggregated patch-level predictions to slide-level predictions using two multiple-instance learning techniques: Max-Pooling (Max) and Dual-Stream Multiple-Instance Learning (DSMIL). Given that optimal hyperparameters for the MIL models may differ based on the representations learned by various SSL methods, we conducted grid searches over learning rates [1e-2, 1e-3, 1e-4] and weight decays [1e-2, 1e-3, 1e-4], yielding 9 hyperparameter combinations per MIL model. This ensured a fair evaluation after optimizing each method’s settings. The MIL models were trained for 50 epochs with a 5-epoch warmup and cosine annealing learning rate schedule. To select the best checkpoint for each representation, we split the training sets into 75% training and 25% validation partitions. The checkpoints with the highest validation set performance were chosen for final evaluation on the held-out test set (Table 2 and 3).

4.3. Ablation Study

To investigate the influence of the core components of our CluSiam model, we conducted an ablation study on the Camelyon16 dataset. We created and trained several model variants with different combinations of the invariance loss term, cluster loss term, and Gumbel noise to evaluate their impacts.

The stop gradient operation is crucial to prevent our model from collapsing, similar to its role in SimSiam. Additionally, incorporating Gumbel noise when assigning clusters is critical. Without this noise, the cluster assigner will collapse early in training, assigning all samples to one cluster. This collapsed state becomes equivalent to SimSiam training, as a single cluster removes the off-diagonal elements required to calculate the cluster loss ${ℒ}_{\text{cluster}}$. The gradient from ${ℒ}_{\text{cluster}}$ thus becomes zero in this situation. However, by introducing randomness in cluster assignment, the noise prevents updates from concentrating solely on the most probable cluster dimensions, thereby preventing early collapse. Together, the stop gradient and noise allow our model to escape these trivial single-cluster solutions, enabling effective joint optimization of ${ℒ}_{\text{inv}}$ and ${ℒ}_{\text{cluster}}$. As evidenced in Table 4, CluSiam outperforms the other models. We observed a collapse in the joint model’s cluster assigner, with one dimension consistently dominating the others. This is likely because only the largest dimension receives updates, continually amplifying its magnitude and leading to a single dominant cluster. The joint model was expected to mirror SimSiam’s performance, as its training dynamics should be identical to SimSiam after the assigner’s collapse. However, its performance was inferior compared to SimSiam. This discrepancy may originate from the cluster assigner offering suboptimal initialization for the subsequent pure SimSiam training before its collapse.

In the aftermath of our ablation studies on the new modules introduced to SimSiam, we delved into investigating the pivotal roles of original SimSiam components in ensuring training stability. A prime area of focus was the scale of inputs to the cluster assigner, as it might significantly influence this stability. The SimSiam projector, which interleaves $ℬ𝒩$ layers between its linear layers and concludes with a $ℬ𝒩$ layer, could be foundational for the clustering module’s effectiveness. To empirically assess the role of controlled input scaling, we devised an experiment on two different projectors. The first projector, originating from SimSiam, concludes with a $ℬ𝒩$ layer. The second projector, from BYOL, ends with a linear layer. We began by replacing the SimSiam-style projector with the BYOL variant in our CluSiam model, leading to the creation of the CluBYOL model with the cluster assigner module integrated into the BYOL architecture. The CluBYOL was initially trained using the BYOL-style projector. Subsequently, we utilized the SimSiam-style projector and undertook another training round for the CluBYOL model. Notably, both situations using the BYOL-style projector resulted in collapse, akin to the joint model depicted in Table 4, with a single cluster predominantly emerging. As highlighted in Table 5, the SimSiam-style projector, characterized by its concluding with a $ℬ𝒩$ layer, is instrumental in preventing such collapses.

To further analyze the behavior of the clustering module, we examined the impact of the exploration space size, denoted as $K$. When Gumbel noise is introduced, some samples exhibit the highest probability of remaining in the cluster with the highest output value, while also possessing a high probability of transitioning to a nearby or similar cluster. Some samples that are difficult to differentiate might be allocated with near-equal probabilities across multiple centroids. A larger value of $K$ results in more refined clustering. Conversely, when $K$ is small, achieving definitive clustering becomes difficult. For instance, with $K$ set to 3, the cluster assigner can allocate some hard-to-distinguish samples into a third cluster, roughly equidistant from the first two. Yet, with $K$ set to 2, the assigner can only place samples into one of the two clusters. Importantly, when $K$ is 1, the model essentially becomes SimSiam since its loss function and backpropagation are equivalent to those in SimSiam. In this scenario, the cluster assigner lacks the flexibility to differentiate between samples by assigning them to different clusters.

In our experiments, we assessed the influence of $K$ on model behavior by training two models with exploration spaces of $K=10$ and $K=100$. Using a Top-1 KNN classifier’s F1 score for patch-level performance evaluation, both models displayed comparable classification performance, as depicted in Figure 4. Despite this similarity in classification, their clustering behaviors were distinct. The model with $K=10$ exhibited a larger number of clusters and greater fluctuations in both cluster counts and the Rand Index. In contrast, the larger exploration space of $K=100$ allowed the assigner to stabilize on more definitive assignments faster. This difference underlines the influence of $K$ in CluSiam’s clustering. Specifically, both the number of clusters and the Rand Index fluctuate more with $K=10$. This limited exploration space prevents highly fine-grained cluster assignments. In contrast, a larger exploration space of $K=100$ provides more granularity for refined clustering actions, enabling the assigner to stabilize on more definitive assignments rapidly. The cluster count also becomes more consistent with $K=100$. These observations highlight the impact of $K$ on CluSiam’s clustering dynamics.