\section{Introduction}

In medicine, large unlabeled image datasets vastly outnumber curated datasets. Furthermore, when labels are available, discrete labels are often employed to describe a continuous spectrum of conditions. There is a lot of interest in developing unsupervised methods to uncover the hidden, or latent, structure of these disorders, in order to identifying disease subtypes or novel classess of disease (see [1] for a review). Unsupervised learning, however, is a challenging task to solve even in domains where large and curated datasets are available. Many diverse approaches have been recently proposed, however they commonly rely on extremely large amounts of data. I norder to render such approaches more practical, a possible solution could be the development of data-efficient unsupervised methods or an increased reliance on self-supervised or unsupervised pretraining followed by supervised fine-tuning, which would drastically reduce the need for labeled samples. In thei context, a number of papers have tackled the problem of clustering using deep learning architectures that find a nonlinear mapping from the input space to the feature space, where clustering is performed (see “Related Work” section).
In this work, we propose an unsupervised deep learning probabilistic approach that projects the inputs to a learned latent space with “prototype” points, such that the low-dimensional representations of the inputs are useful both for their reconstructions and for subsequent tasks such as classifications. Each “prototype” point is associated with a cluster in the latent space. During reconstruction, we condition the decoder on both a) the sample latent embedding and b) a soft cluster assignment, which relies on the similarity of the sample embedding to all prototype vectors. We assume that downstream tasks, such as e.g. classification, are simplified by this conditioning, which provides more context on the relative position of the sample embedding in the latent space. In this context, the main challenge is defining a strategy to learn these prototype vectors, which should act represent subsets that cluster together in the latent space.
The objective of this work is to augment the VAE framework to also learn prototype vectors, hence obtaining a better organization of the latent space. The prototypes are learned using a momentum update during training and are stored in the memory of the model. This conditional information about similarity during decoding should carry all the information needed for reconstructions and implicitly impose a structure in the latent space. We demonstrate that this approach is more effective than two-step approaches that first train a vanilla VAE and then cluster its latent space, for example, by using KMeans. We show that our end-to-end approach outperforms these and other deep clustering baselines on the MNIST handwritten digit dataset  \cite{mnist}as a benchmark and on the Pneumonia dataset from the MedMNIST collection \cite{medmnistv2} as an example in the field of biomedical images. These performance differences are quantified in terms of different metrics, related to both the alignment between the label clusters and the classification obtained using kNN, or linear classification to measure the capacity of the latent space to capture the fundamental characteristics of the encoded samples.

\section{Related Work}
General trends in unsupervised and self-supervised deep tend to be dominated by large networks trained over hundreds of thousands of images. Recent examples include student-teacher approaches such as DINO  \cite{dino}, based on knowledge distillation without the need for labels, and contrastive-based approaches such as \cite{simclr,simclr2,moco}, where models are trained through a contrastive loss and a similarity metric to position similar samples close in the latent space and different ones very far apart and make heavy use of different augmentations to generate positive pairs. However, due to the intrinsic properties of medical image datasets, these methods poorly translate to the medical context. Knowledge distillation approaches such as DINO rely on the use of data-hungry architectures like vision transformers \cite{vit} or ResNet50 \cite{resnet}], whose power can be harnessed only when “big data” is available. Contrastive-based approaches can tackle label scarcity by employing augmentations, but in medical applications, geometrical augmentations such as mirrors, zooms, and color distribution shifts cannot capture the full pathology diversity.
Furthermore, approaches like SimCLR \cite{simclr} need a very large batch size to consider negative sampling, resulting in high computational demands that can only be met with nonstandard hardware. This problem could be partially solved using more efficient methods such as MOCO  \cite{moco}where a momentum encoder is used for negative sampling storing low-dimensional representations of the precedent batches. Even if all these approaches are appealing and promising when tested on large-scale datasets, their direct translation to data from the medical domain is still not straightforward.
Smaller models could be useful when data are scarce and are mainly built on self-supervised pretraining using autoencoders or variational autoencoders. In \cite{deepclustervariationautoenc} ], the authors trained a variational autoencoder using information from a Gaussian mixture model (GMM) employed as an unsupervised clustering algorithm. Including this information in the loss function, they were able to drive the model to learn a well-separated latent space. In  \cite{vade} ], the authors replaced the VAE prior with a GMM-based prior enforcing the same type of separation in a more direct way. In \cite{deepkmeans}.
], a differentiable version of the K-Means algorithm based on the softmax function (instead of argmax) is proposed and embedded in a neural network. All these methods propose different solutions to the clustering problem by exploiting classification with soft labels or even reconstruction of samples, in most of cases using benchmark datasets as proof-of-concept demonstrations. In general, one of the main open issues is how to inject information about cluster labels, which are usually obtained using nondifferentiable functions, into the model. In our work, we combine soft labels obtained by tempered softmax similarity to inform the decoder about clustering and use nondifferentiable functions to update the prototype values.
Other unsupervised approaches extended those and other ideas to the medical field. Examples are the Mixture of Experts (MoE)  \cite{moe}, where a clustering network proposes a cluster and a mixture of experts learn how to reconstruct the data related to single-cell data, or \cite{CNN-Kmeans-medical} where a CNN is combined with a KMeans approach to classifying an image dataset with different modalities. Finally, a translation of contrastive learning, which also uses multi-instance learning, has achieved results comparable to fully supervised state-of-art algorithms \cite{micle}.
Again, however, in many clinical/medical domains, these strategies might be suboptimal because they rely on nondifferentiable methods, and require a large amount of data or domain-specific knowledge necessary to incorporate biologically/medically aware data augmentations.

\section{Materials and methods}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.6\linewidth]{images/overview_architecture_text.png}
    \caption{Overview of VAESim architecture. We use  $\epsilon$ to denote sampling with the reparametrization trick, $sim$ denotes the similarity measure, and $\parallel$ 1- the concatenation operation.}
    label{fig:overview}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.6\linewidth]{images/overview_update_text.png}
    \caption{
    Graphical description of the $Q$ matrix update mechanism. See Algorithm \ref{alg:update} for details. We use the similarity between prototypes and sample embeddings to assign each sample embeddings to assign each sample embedding to the cluster with the most similar prototype vector. Then, we compute the mean of all $z$ embeddings in the same cluster and use the exponential moving average between them to update the prototypes matrix.}
    \label{fig:overview}
\end{figure}

The main goal of VAESim is to augment the vanilla VAE by The main goal of VAESim is to augment the vanilla VAE by conditioning the decoder not only on the latent embedding of a sample but also on its soft cluster assignment in latent space. This involves two main challenges:  learning a descriptive latent space that contains sufficient information to reconstruct the input samples and 2) learning the prototypical vectors of the latent space. The VAESim model is composed of three elements: the encoder  $e$, the decoder $d$, and the prototype matrix $Q^{K\times q}$.

First, we encode the $p$-dimensional vector input $x$ with a neural network $e: \mathbb{R}^p \rightarrow  \mathbb{R}^q\times\mathbb{R}^q $ into a $q$-dimensional latent space $z$. Then, we sample a latent representation from the estimated posterior, modelled as a Gaussian, $z \sim \mathcal{N}(z_{\mu},e^{\frac{z_{\log{\sigma}}}{2}})$.
We then compute the similarity between the latent embedding and the prototypical vectors using cosine similarity, and perform soft cluster assignment by computing the tempered softmax over the similarity measures.

\begin{equation}
    \label{eqn:cluster function}
	c = softmax(\frac{Qz}{\tau})

\end{equation}

We use $\tau$ to denote the temperature parameter and $c$ to denote the cluster assignment. The latent representation $z$ is passed together with the cluster assignment $c$ to the decoder $d$ which maps back $\tilde{x}=d(z,c)$ into the input space to reconstruct the original input sample. We augment the standard ELBO loss of the vanilla VAE with an optional term that enforces orthogonality in the prototype representations.

The prototype vectors are updated following Algorithm \ref{algorithm}. The differentiating steps of our method are the \textit{cluster} equation and the \textit{update} function. The \textit{cluster} equation Eq. \ref{eqn:cluster function} produces soft assignments given the latent embeddings $z$ and the prototypes $Q^{K\times q}$. The \textit{update} function presented in Algorithm \ref{alg:update} is one of the core building blocks of our approach. At each training iteration, we assign all samples in the training dataset to a cluster. Each prototypical vector is updated with an exponential moving average with the mean of all sample embeddings assigned to its respective cluster.

\begin{algorithm}[H]
  \caption{Training VAESim model}
  \label{algorithm}
  \begin{algorithmic}
    \item[\algorithmicinput]{A set of N training images $X^{N\times p}$.}
    \STATE $\tau: 1.$ linearly decreases to 0.01 in half of the epochs.
    \STATE Initialize prototypes matrix $Q$.
    \FOR{$x$ in the training set}
     \STATE optimizer.zero\_grad()
      \STATE $z=e(x)$
      \STATE $c=softmax(\frac{Qz}{\tau})$
      \STATE $\tilde{x}=d(x,c)$
      \STATE loss=$\mathcal{L}_{recon}(x,\tilde{x})+\beta KL_{div}(q_z,p_z)$
      \STATE loss.backward()
      \STATE optim.step()
      \STATE Update(Q,z,c)
    \ENDFOR
  \end{algorithmic}
\end{algorithm}

\begin{algorithm}[H]
  \caption{\textit{Update} function}
  \label{alg:update}
  \begin{algorithmic}
    \item[\algorithmicinput]{$Q$,$z$,$c$}
    \STATE $\eta: 0.95$
    \FOR{$z_i$,$c_i$ in z,c}
     \STATE Sample cluster label $k_i \in [0,K-1 ]$: $k_i \sim multinomial(c_i)$.
    \ENDFOR

    \FOR{$i=0\dots K$}
     \STATE  Update the i-th column of $Q$ matrix using vectors assigned to it
     \STATE  $\bar{z_{i}}$=mean($z_i$) \textbf{if} $k_i==i$
     \STATE $Q_i=\eta \bar{z_{i}} +(1-\eta)Q_i$
    \ENDFOR
  \end{algorithmic}
\end{algorithm}

\subsection{VAE}
The main advantage of a variational autoencoder compared to standard autoencoders is that we impose a well-structured and continuous latent space, allowing us to sample meaningful examples from the latent space. Instead of simply projecting the data points into a manifold and then reconstructing them, in variational autoencoders the encoder neural network tries to approximate the joint probability $p(z|x)$ (which is computationally intractable) through probabilistic encoding $q_z(z|x)$. The encoder outputs the mean and the log-variance of a multidimensional Gaussian distribution. From that distribution, a value is sampled and passed to the decoder which tries to map back from $z$ to $\tilde{x}$.
The model is optimized by enforcing at the same time the quality of reconstruction through an MSE loss and a regularizer term that forces closeness of the Gaussian distribution to a chosen posterior which is a standard normal distribution. The latter term is usually a weighted Kullback-Leibler  (KL) divergence.
Therefore, the model is optimized to reconstruct items that have a latent distribution that resembles as much as possible a standard multivariate Gaussian distribution.
The loss function is therefore $\mathcal{L}=\|x-\tilde{x}\|^2+\beta KL(q_z,\mathcal{N}(0,1))$.
Also, in this work, we modify the decoder network to take as input both the latent representation and a condition vector that is a measure of similarity across prototypes for each sample.

\subsection{Prototype Matrix Q}
When the model is initialized, the required number of clusters $k$ and the dimension of the latent space $p$ are initialized \textit{a priori}. The $Q$ matrix  should be a matrix with shape ($k$, $p$) having $k$ rows, each one with $p$ dimensions.
At the first forward pass, $k$ elements from the first batch are chosen randomly and set as rows of the $Q$ matrix.
In the early phase of training, the values of the prototypes are not meaningful so we use a tempered softmax of the similarity with a temperature schedule. For the first quarter of the training steps, the temperature is set to a high value ( $>1$) flattening the distribution of similarity across different prototypes. While the model’s ability to reconstruct images increases, the temperature decrease, results in harder conditioning vectors.
The key point related to the prototype matrix $Q$ is the way it is updated. Instead of relying on gradient descent, the values of the $Q$ matrix are computed separately from  the computational graph by using a momentum update of all samples belonging to a specific cluster, i.e.s all the samples whose z is most similar to a specific column of the $Q$ matrix.

At each iteration, the batch ``cluster labels'' are sampled from (or obtained by argmaxing) the categorical distribution $c$. These labels are used to compute the \textit{update} function for the $Q$ matrix.

\subsection{Evaluation}
The regularization constraint imposed by the KL divergence in the loss function tends to collapse all inputs as close as possible in the latent space. Thus, using classical metrics for unsupervised clustering could be misleading due to the commonly used normalization by the intercluster distance. Instead, we produce a visually plausible representation in 2D using t-SNE to observe how inputs are arranged in the latent space. Successively, we proposed three approaches for performance evaluation, all based on the downstream classification task. model training is therefore fully self-supervised way, and some of the available labels are used for comparison in the evaluation phase. We apply three standard methods for learning the downstream task: a statistical mapping approach between cluster labels and real labels, the kNN (k-nearest neighbors) algorithm and the training of a linear layer with categorical cross-entropy as the loss function.
In the statistical mapping approach, clusters are defined by argmaxing the similarity between samples and prototypes. Each cluster label is then associated with the most frequent label present in the related cluster, as a a way to quantify how well each prototype captures the key elements of a class. In inference, this mapping is used to predict the final class.
In the second approach, a subset of the training set is used to compute the "memory bank" for the kNN algorithm. In inference, each sample of the test set is compared with each element stored in the memory bank by computing the pairwise  euclidian distance  . The predicted label is defined as the mode of the label of the $k$ closest samples.

In the third approach, a linear layer maps from $z$ to the number of possible classes. This linear layer is trained with the Adam optimizer and a learning rate of $3*10^{-4}$ for 200 epochs over a subset of the training set, with categorical cross-entropy as the loss function.

\subsection{Baselines}
 In order to compare our probabilistic framework to classical two-step approaches, we trained a VAE build to be as similar as possible VAE we use to reconstruct images. Then, we performed a KMeans evaluation to produce cluster labels and train kNN and a linear classifier on the latent space as well as in other experiments. In this way VAE+KMeans can be evaluated on the same metrics used for our approach.   We also evaluated two other approaches, namely VaDE \cite{vade} and GMVAE \cite{gmvae} using default parameters. Each experiment was run $10$ times to obtain mean accuracies as well as standard deviations. It shoul dbe noted that reported baseline accuracies may not overlap with those shown in original papers because of different splits in the training, random seeds and, in particular, different ways to map between cluster labels and ground truth. Specifically, both GMVAE and VaDE rely on a measure of clustering accuracy defined by mapping $N$ cluster labels to $N$ real labels using the Hungarian matching algorithm \cite{hungarian}, which therefore requires the number of clusters to be indentical to the number of distinct labels. We choose to use a different mapping approach between cluster and actual labels to allow more freedom in the choice of cluster numerosity, hence allowing  the potential discovery of finer structures while maintaining the possibility of many-to-one mapping.
As comparison, we also report state-of-art results achieavable with transfer learning from ResNet50 pretrained on ImageNet.

\subsection{Experiments}

We tested our model on two datasets: a classical dataset of handwritten digits (\cite{deng2012mnist}), and a medical benchmarking dataset, PneumoniaMNIST, part of the MedMNIST dataset \cite{medmnist,medmnistv2}, which is a collection of medical imaging datasets arranged for benchmarking.
Pneumonia was chosen as a biomedical example, where chest X-ray images are assigned two possible labels: "healthy" and "pneumonia". This dataset is composed of 5856 images (training test split: 4708/1148). All code is written in PyTorch [20] and trained on a server with an A6000 GPU and 512 GB RAM.

Code for VAEsim implementation and other algorithms used as baselines can be found at \url{https://github.com/matteoferrante/VAESIM.git}

\section{Results}
Our main hypothesis is that our end-to-end deep clustering approach would improve performance in downstream classification tasks as compared to a two-step approach as well as other baselines. For the VAE+KMeans baseline, we train a VAE with architecture as close as possible to the ones used in our approach ( \textit{latent dimension} $32$ and \textit{batch size} $2048$).

Then, we used K-Means to evaluate the optimal number of clusters using the elbow method with the yellow-brick library \cite{yellowbrick}. VaDE and GMVAE were used as baselines for deep clustering. All results are described in `. Our model consistently outperformed (or performed on par in all metrics) within the baselines investigated, with improvements over the two-step approach ranging from $+13\%$ up to $+17\%$ on the MNIST dataset and from $+4\%$ to $+17\%$ on the Pneumonia dataset. The most relevant improvements are on the kNN accuracy measure, suggesting that our modification to the VAE framework drives the model to learn a space better suited for kNN estimations, i.e. positioning similar examples closer to each other in the latent space as compared to a standard VAE+KMeans.

\begin{table}[]

\label{tab:results}
\resizebox{\textwidth}{!}{\begin{tabular}{llllll}
\hline
\multicolumn{1}{|l|}{\textbf{Approach}} & \multicolumn{1}{l|}{\textbf{dataset}} & \multicolumn{1}{l|}{\textbf{statistical acc}} & \multicolumn{1}{l|}{\textbf{knn acc}} & \multicolumn{1}{l|}{\textbf{linear acc}} & \multicolumn{1}{l|}{\textbf{\begin{tabular}[c]{@{}l@{}}resnet50 acc\\ (supervised)\end{tabular}}} \\ \hline
VAE $\pm$ KMeans                              & MNIST                                 & 0.68 $\pm$ 0.01                                     & 0.80 $\pm$ 0.004                            & 0.68 $\pm$ 0.04                                & 0.99                                                                                \\[2pt]
GMVAE                                   & MNIST                                 & 0.76 $\pm$ 0.01                                     & 0.91 $\pm$ 0.001                            & 0.76 $\pm$ 0.01                                & 0.99                                                                                \\[2pt]
VaDE                                    & MNIST                                 & 0.81 $\pm$ 0.05                                     & 0.964 $\pm$ 0.001                           & \textbf{0.83 $\pm$ 0.008}                      & 0.99                                                                                \\[2pt]
\textbf{VAEsim}                           & MNIST                                 & \textbf{0.83 $\pm$ 0.01}                            & \textbf{0.970 $\pm$ 0.001}                  & 0.81 $\pm$ 0.003                               & 0.99                                                                                \\[2pt] \hline
VAE $\pm$ KMeans                              & Pneumonia                             & 0.631 $\pm$ 0.01                                    & 0.613 $\pm$ 0.014                           & 0.451 $\pm$ 0.033                              & 0.85                                                                                \\[2pt]
GMVAE                                   & Pneumonia                             & 0.625 $\pm$ 0.001                                   & 0.627 $\pm$ 0.041                           & 0.429 $\pm$ 0.061                              & 0.85                                                                                \\[2pt]
VaDE                                    & Pneumonia                             & 0.625 $\pm$ 0.001                                   & 0.747 $\pm$ 0.012                           & 0.590 $\pm$ 0.116                              & 0.85                                                                                \\[2pt]
\textbf{VAEsim}                           & Pneumonia                             & \textbf{0.677 $\pm$ 0.001}                          & \textbf{0.778 $\pm$ 0.024}                  & \textbf{0.671 $\pm$ 0.010}                     & 0.85           \\[2pt]
\end{tabular}}

\vspace{2mm}

\caption{
Results on MNIST and PneumoniaMNIST datasets. The results are compared to different baselines and evaluated with accuracy measured with the statistical mapping approach to evaluate overlap between cluster labels and real labels, kNN to measure the closeness between similar samples in the latent space, and the classification done with a linear layer. The last column refers to the accuracy of a fine-tuned model from ResNet50 pretrained on ImageNet.}
\end{table}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.6\linewidth]{images/prototypes_mnist.png}
    \caption{Prototypes for MNIST. These are the reconstruction of the $Q$ matrix rows.}
    \label{fig:mnistprototypes}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{images/tsne_mnist.png}
  \caption{tSNE visualization of the latent space. Left: colors are real labels. Right: colors are cluster labels, with superimposed cluster prototypes.}\label{fig:tsnemnist}
\end{figure}

Examples of the \textit{prototypes} are visible in Fig \ref{fig:mnistprototypes}. Qualitatively, they clearly represent different ways to write digits and this pattern emerges from the update function during training. The  Fig \ref{fig:tsnemnist} shows a tSNE representation of the latent space with both the real labels and the cluster labels.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.7\textwidth]{images/pneumonia_prototypes.png}
    \caption{Prototypes for Pneumonia. These are the reconstruction of the $Q$ matrix rows.}
    \label{fig:pneumoniaprototypes}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.7\textwidth]{images/tsne_pneumonia.png}
  \caption{tSNE visualization of the latent space. Left: colors are real labels. Right: colors are cluster labels, with superimposed cluster prototypes.}\label{fig:tsnepneumonia}
\end{figure}

Prototypes and t-SNE representations for Pneumonia datasets are visible in the Fig \ref{fig:tsnepneumonia} and Fig \ref{fig:pneumoniaprototypes}.
The prototype representations include both some artefactual structures and some realistic X-ray chest images, where in many images pneumonia-related opacities are visible.

\section{Discussion}

In this work, we proposed a probabilistic framework that modifies the standard training of a variational autoencoder to generate representations of prototypes or modes present in the dataset. After encoding, the representations are compared with a memory bank of those prototypes, and a measure of similarity is passed to a conditional decoder, together with the latent representation, to generate reconstructions. The way we update the values of the prototypes is detached from the computational graph, allowing us to aggregate and update the information using nondifferentiable functions. As proof of concept, we choose an exponential moving average of the closest samples to the cluster prototype. At each step, the new value is computed from the actual value and average of all the samples in batch size for which this prototype is at a minimum distance.
In our approach, in the first part of the training, the focus of the encoder is to indirectly generate robust prototypes. It is an indirect process because the update of Q relies only on how close and similar same-class labels are in the latent space. Forcing the decoder to deal with the conditioning vector and using this updating procedure, we enforce an organization in the latent space that is better than the ones resulting from standard variational autoencoders. During the training process, the temperature  of the similarity conditioning vector decreases, resulting in harder assignments. When the prototypes became more stable, the encoder focuses more on exploiting features in terms of differences between samples and prototypes. We showed that our continuous conditional learning with decreasing temperature improves the latent space organization with representation, which is also useful for the downstream classification task and results in better performances as compared to two-step approaches. Also, our approach outperforms the baselines of other approaches based on deep clustering variational autoencoders. All other VAE proprieties are unchanged, hence allowing to sample from the latent space and reject bad samples based on their similarity with a specific cluster. For example, looking at cluster prototypes obtained for the model trained on MNIST, one could select clusters that represent a digit or a specific way to write a digit and sample from the latent space random values, possibly keeping only the ones that surpass a specific threshold of similarity with that cluster, to reconstruct similar digits. In the same way, if able to interpret prototypes for X-ray imaging or other medical imaging, one could reconstruct specific modes of healthy or pathological images to investigate their proprieties. Importantly, our model is not targeted to maximum accuracy because it is trained for a different objective (image reconstruction). Still, as an emergent property, we obtain a well-organized latent space that can also be exploited for a downstream task. Our model can perform simultaneous reconstruction and clustering. This could be useful, for example, in guided data augmentation or in patient stratification beyond what is visible to the naked eye of the clinician.

It should be noted that several hyperparameters need initialization. In particular, the most challenging one is the number of prototypes. In  \ref{sec:appendix} An investigation of the hyperparameter effect is presented, where we acknowledge that the computational efforts necessary to find the best hyperparameters for any given task could represent one of the main limitations of this model. Imposing a high number of prototypes may improve performance, however, each cluster contains a sufficient number of samples, hence requiring a larger batch size. Finding the optimal way to initialize the Q matrix in terms of values and dimensions will be the object of future investigations.

\section{Conclusion}
 In this work, we proposed a probabilistic framework that improves latent space generalizability for reconstruction and downstream tasks by learning how to reconstruct images through latent representation coupled with a measure of similarity with prototypes that represent the modes present in the dataset. We dynamically update representations during the training, forcing the model to implicitly take into account features that enhance similarities and differences between samples and prototypes. We demonstrate that our approach is better than a two-step approach basedon the combination of a VAE and KMeans, and training these models in an end-to-end fashion could lead to up $+10\%$ in performance on average. We also outperformed baselines based on different modifications of the VAE framework for deep clustering on most of the metrics that we defined. Also, our model produces a latent space that generalizes well when used for the downstream task. We evaluated performance by examining the ability of cluster labels to recall specific classes, the organization in latent space probed by the use kNN algorithm for classification and its separability using a linear classifier. We also reached performances similar to the ones obtained fine-tuning very large models in a supervised way by starting with a self-supervised pretraining of our modified VAE.

\bibliographystyle{splncs04}


\section{Introduction}

\subsection{Motivation}
So let's say we're trying to develop an autonomously driving car as many auto companies are trying to do now. There are many different problems to consider if we're to give this "driver" a low chance of ever causing problems that exceed an acceptable degree of damage. We could try to train this driver based on what human drivers do in certain scenarios. This method could be employed through the use of a neural network. The process would entail feeding information about the environment to the driving agent and regressing its decisions against many different human decisions in the same environment. These methods work well for scenarios that are well covered by the data, but let's say this robot driver were to encounter a scenario that the data did not cover well if at all. For example many drivers stop for bikes or pedestrians going across a crosswalk, but what about people walking their bikes? This scenario happens less frequently and if our data set or environment did not adequately detail the actions associated with seeing a bike-walker going across a crosswalk, the driving agent will still make a confident decision that could be correct, or wrong. If the robot driver decides the bike-walker is nothing more than debris, it will not stop for them. This type of agent would not be very desirable for this problem since across many instances of autonomous vehicles in uncovered situations by the data, we could expect many accidents with irreparable damage. Thus a new model may be more appropriate.

Instead of trusting the data to account for corner cases in the task you are training for, we should manipulate the data to create the worst possible responses for the neural net we have trained. Once a net is trained, if a noise pattern is applied to whatever input in such a way that it minimizes performance within certain limitations for a given noise additive, then the neural model will typically fail miserably. In image recognition for example, if a neural net is used to categorize the contents of an image, it will perform very well on organic data, which has no manipulations. Once the images are doctored to cause the net to fail, it will select an incorrect category with high confidence when a human would easily identify the image correctly. What's useful about this process is that it identifies weaknesses in the model which would otherwise stand if the net continued to regress against regular labeled data, but by regressing against the doctored cases alongside regular data, the neural net will reach a much more robust solution to the task. This type of training is well understood in the machine learning field, but there is another type of training  that has unique benefits which use these methods and do not have well-understood consequences.

Going back to the self-driving car example, instead of using a data set based on what people experience in everyday driving, let's try and simulate the driving experience on a computer where the consequences for failure are minimal. In this simulation, we could give the robot driver incentives and penalties based on how their performance relates to our goals. Now we are in the realm of reinforcement learning. In this scenario, we don't need data from other drivers, we aren't assuming that we have optimal decisions associated with given scenarios, and there is more choice as to what the agent encounters and how it changes its behavior as a result. Unfortunately given our limited perception of all the potential failures for this robot driver, we cannot manually attend to its goals and expect it to reach a comprehensive model of a good driver. Not only would explicitly setting the goals like where to go and how to get there for the driving agent be tedious and time-consuming but our own bias as to what a driver would need to understand would flaw the overall results of the training. Instead, we could assign another agent the task of setting these training goals with its own goal of causing as much damage and taking as much time as possible, the opposite of our driving agent's goals. In this scenario, if the driving agent had a glaring weakness for driving in certain situations, such as bike walkers crossing the road, the goal-setting agent would exploit them, asking the driving agent to come up to as many crosswalks with bike walkers as it could until the driving agent learned to drive well in these scenarios. Once the driving agent corrected its driving policy the goal-setting agent would lose interest in those road conditions searching for new ones to create problematic results. This process would continue to force the driving agent to patch the holes in its policy, and if it eventually reached a stalemate or a saddle point with the goal-setting agent in terms of performance, then the resulting driving policy would perform very reliably in actual road settings. A major question that comes up, is if these two opposing systems will reach such a point. This depends upon the algorithms behind each of the learner systems and the task they optimize around. For a particular type of algorithm, and problem this convergence is guaranteed~\cite{Chambolle2011,Chambolle2016} and we hope to observe this as we implement these models in simplistic scenarios.

\subsection{Problem Description}

There are several key characteristics of the reinforcement learning design or Markovian Problem setting. The first is that the task and model interface is modeled as a cyclical relationship between an environment and an agent (the algorithm that is being trained). The environment is a set of states that the agent responds to by selecting an action from the set of potential actions. Thus the cycle goes as follows: There is some initial state which is taken by the agent and in return, the agent selects an action. This action is fed into the environment, and the environment transitions the agent into a new state. The agent then selects a new action based on the new state, and the environment changes the state again. This process can continue for some set length or indefinitely. Along with a new state, the environment also returns a reward as a function of the state that the agent was in alongside the action it took in response. This reward is what the agent optimizes around and there are a variety of different methods for doing this.

One of the most successful reinforcement learning methods is ``Policy Gradient'' which optimizes an equation that represents how the reward received from the environment changes as a function of the agent's decision-making algorithm. In many cases, the agent's decision-making is characterized by a set of parameters that correspond to the preferences it has towards certain actions given a set of features in the state. Once an equation that models this relationship is derived, the algorithm's preferences can be optimized around the average reward earned in the environment. This process entails adjusting the agent's parameters proportionally to how they change the reward, which is approximated by the change modeling function.

When two policy gradient algorithms have opposite reward functions for the same environment, they are in an adversarial setup where each will optimize their strategies to directly conflict with the other algorithm's policy. If the two policies reach a balancing point where either one of the agents changing their policy will cause them to perform worse, they are said to be at a \emph{saddle point} or an equilibrium~\cite{Behzadian2021,Ho2021a}. What we're interested in is whether or not reaching this point is guaranteed based on the algorithms and the problem. In particular, we would like to know if certain policy gradient algorithms converge in a variety of complex scenarios which are modified versions of simplistic ones.

\section{Static Problem}
By viewing the relationship between the environment and the agent as a game where the agent maximizes reward, and the environment minimizes it, where the agent is a probability distribution of actions and the environment is a probability distribution of states, with a set reward defined for each state action pair, a simplistic Markovian problem is created. Here the probability of encountering different states is not dependent upon actions taken, nor is the probability of taking different actions dependent on the state. This formulation has the following consequences.

The main objective is to solve the following max-min problem for some reward matrix $A \in \mathbb{R}^{n\times n}$ with opposing player vectors $x$ and $y$ of length $n$:
\[
    \max_{x \in \Delta^n} \, \min_{y \in \Delta^n} \; F(x,y)
\]
where
\begin{align*}
    F(x,y) &\;=\; \ x^T A y \\
    \Delta^m &\;=\; \left\{ v \in \mathbb{R}^m \mid \sum_{i=0}^m v_i = 1,\; v \geq 0\right\}
\end{align*}

In order to guarantee that the updated vectors stay in $\Delta^n$, the update formulas take the form:

\begin{align*}
   x^{k+1} &= x^k + \alpha(\bar{x} - x^k)
\end{align*}

\begin{align*}
   y^{k+1} &= y^k + \alpha(\bar{y} - y^k)
\end{align*}

with $x^k, y^k \in \Delta^n$ where $\bar{x}$ and $\bar{y}$ are to solutions to a gradient descent subproblem which which optimizes $x \in \Delta^n$ for an objective function $G$, and $\alpha \in [0,1]$ as a stepsize parameter.

A couple of notes for the optimization methods below. By simply taking vectors $\hat{x}$ and $\hat{y}$ as the mean of the iterates $x_k$ and $y_k$ respectively, a stronger solution to original problem is obtained for sufficiently large numbers of iterations.
Since $F(x,y)$ is well defined, computing its gradients is simple:

\begin{eqnarray}
\nabla_x F &=& Ay \\
\nabla_y F &=& A^Tx
\end{eqnarray}

Also a key condition for a constrained minimum $x \in \Delta^m$ of an objective function $G$ is

\begin{equation}
\forall i,j \quad \left(x_i > 0 \implies \frac{\partial G(x_i)}{\partial x_i} \leq \frac{\partial G(x_j)}{\partial x_j}\right)
\end{equation}
Essentially this means that all the non-zero positions of the optimal solution $x$ must have partial derivatives of G equal to the minimum partial derivative of G in the vector. If this were not the case, then a more optimal solution could be obtained by simply decreasing one of the non-zero positions, and increasing a position which has a smaller partial derivative by the same amount. This will be useful for deriving the update formulas of each method.

\section{Conditional Gradient}

The first method, conditional gradient descent also known as Frank-Wolfe, is characterized by the subproblems:
\[
\bar{x} = \min_{x \in \Delta^n} \left[ -\nabla F(x^k)^T(x-x^k)\right]
\]
\[
\bar{y} = \min_{y \in \Delta^n} \left[ \nabla F(y^k)^T(y-y^k)\right]
\]

Let $G$ be the objective function contained in the minimization for x, then an optimal solution must satisfy (3). Since $\frac{\partial G(x_i)}{\partial x_i} = -\nabla F(x^k)$, only coordinates with maximal corresponding gradient terms can be non-zero. There can be multiple solutions for $\bar{x}$, but by simply committing to one term which has a maximal gradient term, a valid solution can consistently be obtained. Formally that is:

\begin{align*}
   \bar{x} &= \left\{ x \in \mathbb{R}^n \mid  \begin{cases} x_i = 1, \text{ if } i = j \\
                                                         x_i = 0, \text{ otherwise } \end{cases} \right\}
\end{align*}
where
\begin{align*}
    j = \arg\;\max_{i=1,...,n} (Ay)_i
\end{align*}
and the same can be done for y with a slight modification. Since $\frac{\partial G(y_i)}{\partial y_i} = \nabla F(y^k)$, only coordinates with minimal corresponding gradient terms can be non-zero. This makes
\begin{align*}
   \bar{y} &= \left\{ y \in \mathbb{R}^n \mid  \begin{cases} y_i = 1, \text{ if } i = j \\
                                                         y_i = 0, \text{ otherwise } \end{cases} \right\}
\end{align*}
where
\begin{align*}
    j = \arg\;\min_{i=1,...,n} (A^Tx)_i
\end{align*}

Below is a figure of
\[F(\hat{x},\hat{y}) = \text{Reward} \]
\[\max_{x \in \Delta^n} F(x,\hat{y}) =  \text{Upperbound} \]
\[\min_{y \in \Delta^n} F(\hat{x},y) = \text{Lowerbound} \]
given a 3x3 grid representing x and y playing rock paper scissors where a win for $x$ is $+1$, a loss is $-1$, and a draw is $0$, plotted over 1000 iterations using Conditional Gradient Descent, setting $\alpha = .01$.

\includegraphics{figures/conditional_perf.png}

Here is another figure which represents $x$ and $y$'s path through the 3-dimensional probability simplex for Conditional Gradient Descent where $x$ is green and $y$ is red.

\includegraphics{figures/conditional_path.png}

\section{Euclidean Penalized Gradient}

If instead of just minimizing the dot product of the cost gradient and change in position at each step, we solve
\[
\bar{x} = \min_{x \in \Delta^n} \left[ -\nabla F(x^k)^T(x-x^k) + \beta \cdot \| x - x^k \|_2^2\right]
\]
\[
\bar{y} = \min_{y \in \Delta^n} \left[ \nabla F(y^k)^T(y-y^k) + \beta \cdot \| y - y^k \|_2^2\right]~,
\]
where $\beta \in \mathbb{R}_+$ and $G$ is the objective function minimized by $x$ again, then in order to compute an $\bar{x}$ which satisfies (3), $\frac{\partial G(x_i)}{\partial x_i}$ must be found. That is:
\begin{align*}
    \frac{\partial G(x_i)}{\partial x_i} =& \frac{\partial \left( -\nabla F(x_i^k)^T(x_i-x_i^k) + \beta \cdot \|x_i-x_i^k\|_2^2 \right)}{\partial x_i}
\end{align*}
\begin{align*}
      =& -\nabla F(x_i^k) + 2\beta (x_i - x_i^k)
\end{align*}
\begin{align*}
      =& -\nabla F(x_i^k) - 2\beta x_i^k + 2\beta x_i
\end{align*}
Since the optimization condition above needs all non-zero coordinates to have minimal and equal partial derivatives of $G$, and all values of $x$ must be positive, it can be seen that the terms containing $x_i^k$ will determine the values of $i$ for which $x_i$ will be zero. Specifically, if
\[
d_i = -\nabla F(x_i^k) - 2\beta x_i^k
\]
then
\[
(x_i = 0)  \implies \left( d_j \geq d_i \implies x_j = 0\right)
\]
This is true since if there were an $x_j$ greater than zero for which $d_j$ was greater than $d_i$ where $x_i$ was zero, then the partial derivative of $x_j$ would be larger than $x_i$ and its value would be non-zero, which violates the optimality conditions from before. This idea lends itself to a certain method for computing $\bar{x}$. Since we know $\bar{x}$ will be in the simplex we can solve the non-zero $x_i$ to make all their $\frac{\partial G(x_i)}{\partial x_i}$ equal to some $\lambda$ where $\lambda$ is chosen such that the sum of the $x_i$ is equal to one. If $x_i$ is some nonzero coordinate of the optimal solution $\bar{x}$, then it can be solved in the following fashion:
\begin{align*}
    \lambda &= -\nabla F(x_i^k) - 2\beta x_i^k + 2\beta x_i
\end{align*}
\begin{align*}
    \lambda + 2\beta x_i^k + \nabla F(x_i^k) &= 2\beta x_i
\end{align*}
\begin{align*}
    \frac{\lambda + 2\beta x_i^k + \nabla F(x_i^k)}{2\beta} &= x_i~.
\end{align*}
Now since
\[
\sum_i x_i = 1
\]
it follows that
\begin{align*}
\sum_i \frac{\lambda + 2\beta x_i^k + \nabla F(x_i^k)}{2\beta} &= 1
\end{align*}
\begin{align*}
\sum_i \lambda + 2\beta x_i^k + \nabla F(x_i^k) &= 2\beta
\end{align*}
\begin{align*}
n\lambda &= 2\beta - \sum_i 2\beta x_i^k + \nabla F(x_i^k)
\end{align*}
\begin{align*}
\lambda &= \frac{2\beta - \sum_i 2\beta x_i^k + \nabla F(x_i^k)}{n}~.
\end{align*}
By substituting in $\lambda$
\begin{align*}
x_i &= \frac{\frac{2\beta - \sum_i 2\beta x_i^k + \nabla F(x_i^k)}{n} + 2\beta x_i^k + \nabla F(x_i^k)}{2\beta}
\end{align*}
\begin{align*}
 &= \frac{2\beta - \sum_i [2\beta x_i^k + \nabla F(x_i^k)] + n(2\beta x_i^k + \nabla F(x_i^k))}{n2\beta}
\end{align*}
which is $x_i$ in terms of its previous vector $x^k$, $\beta$, and $n$(the number of nonzero coordinates in $x$). Now since it's not known which coordinates in $\bar{x}$ are zero, an easy way to find out is to start by assuming that all the positions in $\bar{x}$ are non-zero and then solve for each $x_i$. If there are negative values in the solution then at least one of the values in the true solution of $\bar{x}$ is zero. To test if there is exactly one zero in the true solution, all that has to be done is to re-solve $x$ excluding the position corresponding to the greatest $d_i$, since that position is the only position that could be zero with only one $x_i$ equal to zero without violating (3). This process of eliminating the positions of $x$ which correspond to the next largest $d_i$ and re-solving can be continued until none of the positions in the solution are negative. In that case the solution is equal to $\bar{x}$. This process can be used to solve $\bar{y}$ as well if its applied to the y-minimization instead.
\[\]
Below is a figure of
\[F(\hat{x},\hat{y}) = \text{Reward} \]
\[\max_{x \in \Delta^n} F(x,\hat{y}) =  \text{Upperbound} \]
\[\min_{y \in \Delta^n} F(\hat{x},y) = \text{Lowerbound} \]
given a 3x3 grid representing x and y playing rock paper scissors where a win for $x$ is $+1$, a loss is $-1$, and a draw is $0$, plotted over $1000$ iterations using Euclidean Penalized Gradient Descent, setting $\alpha = .01$ and $\beta = .5$.

\includegraphics{figures/euclidean_perf.png}

Here is another figure which represents $x$ and $y$'s path through the 3-dimensional probability simplex for Euclidean Penalized Gradient Descent where $x$ is green and $y$ is red.

\includegraphics{figures/euclidean_path.png}

\section{KL Divergence Penalized Gradient}

If $\beta \cdot \| x - x^k \|_2^2$ was substituted with $\beta \cdot \sum_i x_i \log(\frac{x_i}{x_i^k})$ instead, the new minimization would penalize based on KL divergence of $x$ from $x^k$ instead of the Euclidean distance. Thus the new problems to be solved at each step would be:

\begin{align*}
    \bar{x} = \min_{x \in \Delta^n} \left[ -\nabla F(x^k)^T(x-x^k) + \beta \cdot \sum_i x_i \log \left(\frac{x_i}{x_i^k}\right) \right]~.
\end{align*}
\begin{align*}
    \bar{y} = \min_{y \in \Delta^n} \left[ \nabla F(y^k)^T(y-y^k) + \beta \cdot \sum_i y_i \log \left(\frac{y_i}{y_i^k}\right) \right]~.
\end{align*}

If the same approach for finding the solution $\bar{x}$ is used from earlier, where each partial derivative of the minimization function $G$ is set to some constant $\lambda$ then the values of $x_i$ can be derived as follows:

\begin{align*}
    \lambda &= -\nabla F(x_i^k) + \beta (1 + \log\left(\frac{x_i}{x_i^k}\right))
\end{align*}
\begin{align*}
    \lambda &= -\nabla F(x_i^k) + \beta - \beta\ \log(x_i^k) + \beta\ \log(x_i)
\end{align*}
\begin{align*}
    \lambda + \nabla F(x_i^k) - \beta + \beta\ \log(x_i^k) &= \beta\ \log(x_i)
\end{align*}
\begin{align*}
   \frac{\lambda + \nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) &= \log(x_i)
\end{align*}
\begin{align*}
   x_i &= \exp\left( \frac{\lambda + \nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)
\end{align*}
and by solving $\lambda$ such that
\[
\sum_i x_i = 1
\]
\begin{align*}
    \sum_i \exp\left( \frac{\lambda + \nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right) &= 1
\end{align*}
\begin{align*}
    \exp(\frac{\lambda}{\beta})\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right) &= 1
\end{align*}
\begin{align*}
    \exp(\frac{\lambda}{\beta}) &= \frac{1}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)}
\end{align*}
\begin{align*}
    \frac{\lambda}{\beta} &= \log\left( \frac{1}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)} \right)
\end{align*}
\begin{align*}
    \lambda &= \beta\ \log\left( \frac{1}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)} \right)
\end{align*}

By plugging in the solved $\lambda$
\begin{align*}
   x_i &= \exp\left( \frac{\beta\ \log\left( \frac{1}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)} \right) + \nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)
\end{align*}
\begin{align*}
   x_i &= \exp\left( \log\left( \frac{1}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)} \right)+\frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)
\end{align*}
\begin{align*}
   x_i &= \frac{\exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)}{\sum_i \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 + \log(x_i^k) \right)}
\end{align*}
\begin{align*}
   x_i &= \frac{x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 \right)}{\sum_i x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta} - 1 \right)}
\end{align*}
\begin{align*}
   x_i &= \frac{e^{-1} x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta} \right)}{e^{-1}\sum_i x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta}\right)}
\end{align*}
\begin{align*}
   x_i &= \frac{x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta} \right)}{\sum_i x_i^k \exp\left( \frac{\nabla F(x_i^k)}{\beta}\right)}
\end{align*}

Unlike the Euclidean penalized minimization, by simply setting all the partials of $G$ equal to one another, the values of the resulting $x$ are guaranteed to be positive and sum to one. Therefore, there is no need to check how many zeros will be in the true solution, since the first $x$ yielded from the above formula will be in the n-dimensional simplex and meet (3), and thus equal $\bar{x}$. The same procedure can be used for solving $\bar{y}$ by replacing the components from the x-minimization with the components of the y-minimization.
\[\]
Below is a figure of
\[F(\hat{x},\hat{y}) = \text{Reward} \]
\[\max_{x \in \Delta^n} F(x,\hat{y}) =  \text{Upperbound} \]
\[\min_{y \in \Delta^n} F(\hat{x},y) = \text{Lowerbound} \]
given a 3x3 grid representing x and y playing rock paper scissors where a win for $x$ is $+1$, a loss is $-1$, and a draw is $0$, plotted over $1000$ iterations using KL Divergence Penalized Gradient Descent, setting $\alpha = .01$ and $\beta = .25$.

\includegraphics{figures/softmax_perf.png}

Here is another figure which represents $x$ and $y$'s path through the 3-dimensional probability simplex for KL Divergence Penalized Gradient Descent where $x$ is green and $y$ is red.

\includegraphics{figures/softmax_path.png}

\section{Future Work}
Although each of these methods solves the max-min of reward for static distributions of actions and states, this scenario does not adequately encapsulate most Markovian problems. Usually the probability of encountering different states is dependent on the actions taken, and the actions taken should depend on the state the agent is in. There are more complex models known as Distributionally Robust Markov Decision Processes which account for these considerations. Instead of assuming a static distribution of states, a probability distribution of next states is assigned to each state-action pair, and instead of the adversarial environment selecting a distribution of states, it selects a distribution of transition rates across the different state action pairs. The agent also has separate probabilities of selecting different actions in different states. Each of these additions to the problem adds another layer of complexity onto the previously solved static problem, but the iterative methods for solving these problems remain similar.

In order to understand the consequences for making the model more applicable, we would need to spend time looking at the current literature describing these systems, implement the algorithms to solve them, and compare their results against our own theoretical results and see how much they differ. If these general purpose algorithms could have better understood consequences, it would have large implications for the many fields that depend on these models.

\bibliographystyle{plain}


\section{Introduction}
\section{The $W$ boson mass shift and the weak SU(2) gauge coupling constant}
The CDF collaboration recently reported an updated precise measurement of the W-boson mass~\cite{CDF:2022hxs}, which is significantly larger than the Standard Model (SM) prediction and
the world average without this new CDF result~\cite{Workman:2022ynf}.
With this new CDF result, the experimental average of the W-boson mass becomes~\cite{deBlas:2022hdk}
\begin{eqnarray}
(M_W)_{\rm exp} =(80.4133 \pm 0.0080)\, {\rm GeV}, \label{eq:ave_exp}
\end{eqnarray}
which deviates from the SM prediction~\cite{deBlas:2022hdk} by 6.5 $\sigma$:
\begin{eqnarray}
(M_W)_{\rm SM} = (80.3499 \pm 0.0056)\, {\rm GeV}. \label{eq:sm_pre}
\end{eqnarray}

The discrepancy, $\delta M_W \equiv (M_W)_{\rm exp} - (M_W)_{\rm SM} \sim 60$\,MeV
is easily explained by the small vacuum expectation value (VEV) of a $SU(2)_L$ triplet Higgs boson with zero hypercharge~\cite{Strumia:2022qkt, DiLuzio:2022xns, Athron:2022isz, FileviezPerez:2022lxp, Evans:2022dgq}, which only contributes to the $W^\pm$ mass at the tree-level~\cite{Ross:1975fq, Gunion:1989ci, Blank:1997qa, Forshaw:2003kh, Chen:2006pb, Chankowski:2006hs, Chivukula:2007koj}. The mass shift of the W-boson may suggest the existence of the triplet Higgs boson at the TeV scale.

Let us discuss the W-boson mass shift, $\delta M_W \sim 60\,{\rm MeV}$, with the complex $SU(2)_L$ triplet Higgs boson, $\Sigma$, added to the SM Higgs potential. We consider the following potential:
\begin{eqnarray}
 V &\ni&  - m_H^2 |H|^2 + \lambda_H |H|^4 \nonumber \\
  &+& 2  m_{\Sigma}^2 {\rm Tr } \Sigma^\dag \Sigma + (A_{H} H^\dag \Sigma H + h.c.), \label{eq:potential_triplet}
\end{eqnarray}
where we write
\begin{eqnarray}
	\Sigma = \Sigma^a T^a = \frac{1}{2}
	\left(
	\begin{array}{cc}
		\Sigma^3	&  \sqrt{2} X_1^+ \\
		\sqrt{2} X_2^- & -\Sigma^3
	\end{array}
	\right).
\end{eqnarray}
Here, $T^a$ is a generator of $SU(2)_L$, which satisfies ${\rm Tr}(T^a T^b) = 1/2 \delta^{ab}$. By minimizing the potential in Eq.~\eqref{eq:potential_triplet}, the vacuum expectation values for the doublet and triplet Higgs are
\begin{eqnarray}
v^2 =(\mu_H^2 + A_{H} v_T)/(2\lambda_H), \ \
v_T =  \frac{A_{H} v^2}{2\mu_3^2},
\end{eqnarray}
where
\begin{eqnarray}
\left< H \right> = (0, v)^T, \
\left< \Sigma \right> = \frac{1}{2}
\left(
\begin{array}{cc}
 v_T	&  0 \\
 0	& -v_T
\end{array}
\right)
\end{eqnarray}

At the tree-level, the W-boson mass with the triplet Higgs is expressed as (see e.g.,~\cite{Diessner:2019ebm})
\begin{eqnarray}
M_W &\approx&  (M_W)_{\rm SM} \left[ 1 + \frac{1}{2} \frac{\cos^2 \theta_W}{\cos^2 \theta_W-\sin^2 \theta_W} (\alpha T_{\rm tree})  \right], \nonumber \\
\end{eqnarray}
where $\theta_W$ is the weak mixing angle, $\alpha$ is the fine structure constant and $T_{\rm tree}$ is the $T$-parameter at the tree-level:
\begin{eqnarray}
T_{\rm tree} = \frac{4 v_T^2}{v^2} \frac{1}{\alpha}.
\end{eqnarray}
We see that $v_T \sim 3$\,GeV explains the discrepancy, $\delta M_W$.
The loop effects to the oblique parameters, $S$ and $T$, are suppressed as $\sim k^2/16\pi^2 (M_W^2/m_{\Sigma}^2)$ ($k$ is a relevant coupling), which are negligible for $m_{\Sigma} \geq 1\,{\rm TeV}$. For the $U$ parameter, since it comes from a dimension eight operator, the effect is even more suppressed. So far, only $T \simeq T_{\rm tree}$ is non-negligible. This is consistent with the electroweak fit results in Refs~\cite{Strumia:2022qkt, deBlas:2022hdk} (see also \cite{Lu:2022bgw}).

\begin{figure}[htp]
	\centering

	\includegraphics[scale=0.52]{unification_1a.pdf}
	\includegraphics[scale=0.52]{unification_1b.pdf}
	\caption{Renormalization group evolution of gauge coupling constants in SM (dotted), SM + one real triplet Higgs (dashed) and SM + one complex triplet Higgs (solid). The blue, green and orange lines show $U(1)_Y$, $SU(2)_L$ and $SU(3)_C$ couplings, respectively.
	The mass of the triplet Higgs is taken as 1\,TeV}
	\label{fig:1}
\end{figure}

The mixing between the SM Higgs and the triplet Higgs comes from the scalar trilinear coupling in Eq.\eqref{eq:potential_triplet}. The trilinear coupling $A_H$ is related to $m_\Sigma$ as
\begin{eqnarray}
\frac{A_H}{m_\Sigma} = \frac{2 v_T m_\Sigma}{v^2} \approx 0.3
\left(
\frac{m_\Sigma}{1500\,{\rm GeV}}
\right)
\left(
\frac{v_T}{3\,{\rm GeV}}
\right). \label{eq:ratio_ah}
\end{eqnarray}
The ratio, $A_H/m_{\Sigma}$, increases linearly with $m_{\Sigma}$. This relation is almost unchanged even if quartic couplings involving $\Sigma$ are included (unless $m_\Sigma^2 < 0$).
The mixing angle between the SM-like CP-even Higgs and the heavy CP-even Higgs is given by~\cite{Evans:2022dgq}
\begin{eqnarray}
\tan 2 \theta_{\rm mix} &=& \frac{2 A_{H} v}{m_{\Sigma}^2 - 4 \lambda_H v^2} \nonumber \\
&\approx&
\frac{2 A_{H}}{m_\Sigma} \frac{v}{m_\Sigma}
\approx
0.07
\left(
\frac{v_T}{3\,{\rm GeV}}
\right).
\end{eqnarray}
The mixing angle is almost fixed by the VEV of the triplet as long as $m_{\Sigma} \gg v$.
Therefore, the SM-like CP-even Higgs couplings are modified by a fixed amount (up to $\sim\,4\%$) compared to those of the SM Higgs in the heavy triplet region, which is a kind of general prediction.

On the other hand, when $A_H/m_\Sigma \gtrsim 1$, it induces a deeper global minimum other than the electroweak symmetry breaking minimum even if we include quartic couplings containing $\Sigma$.\footnote{For instance, one can consider the direction $\left<h\right>=\left<\sigma\right>=x\  (h \in H, \sigma \in \Sigma^3)$ in the scalar potential. Then, the potential can be written as $V_{\rm eff} \simeq m_{\Sigma}^2 x^2 - 2 A_{3 H} x^3 + l_{\rm eff} x^4$, which has the deeper global minimum for $A_{3 H}^2/l_{\rm eff} > (8/9) m_{\Sigma}^2$.} We expect the vacuum stability puts the upper-bound on $m_{\Sigma}$, which is not much larger than 5\,TeV.

\vspace{20pt}
The presence of the complex $SU(2)_L$ triplet with the zero hypercharge at the TeV scale changes the renormalization group evolution of the gauge coupling constants. In particular, this triplet Higgs only affects the renormalization group equation (RGE) for the $SU(2)_L$ gauge coupling at the one-loop level:
\begin{eqnarray}
\frac{d \alpha_2^{-1}}{d \ln Q_r} = -\frac{1}{2\pi} \left(
b_{\rm SM} + \frac{2}{3}
\right), \label{eq:a2beta}
\end{eqnarray}
where $Q_r$ is the renormalization scale and $b_{\rm SM}=-19/6$.
In Fig.~\ref{fig:1}, the RG evolution of the gauge coupling constants at the two-loop level are shown. We use {\tt PyR@TE 3}~\cite{Sartore:2020gou} to obtain the two-loop RGEs.
Interestingly, as shown in Ref.~\cite{Evans:2022dgq} and Fig.~\ref{fig:1}, the complex $SU(2)_L$ triplet with the mass of $\sim 1$\,TeV leads to rather precise gauge coupling unification. This suggests us to consider a grand unified theory (GUT) at the scale around $10^{14}$\,GeV. It has been shown that $SU(5)$ GUT can successfully be constructed without inducing rapid proton decays~\cite{Evans:2022dgq}.

\section{The quintessence electroweak axion dark energy }

Now what if there were no GUT? Surprisingly, in this case, the observed dark energy may be explained by the potential of the electroweak (EW) axion \cite{Fukugita:1994xx}.\footnote{A proposal to identify the quintessence for the dark energy with a pseudo Nambu–Goldstone boson (axion) is first discussed in \cite{Fukugita:1994xx}, \cite{Frieman:1995pm}, \cite{Carroll:1998zi} and \cite{Choi:1999xn}.}
It was pointed out that the $SU(2)_L$ instantons give a correct magnitude of the observed dark energy in the supersymmetric (SUSY) standard model~\cite{Nomura:2000yk}. Here, we assume that a massless Nambu-Goldstone boson (i.e. EW axion), $A$, couples to the $SU(2)_L$ gauge field as
\begin{equation}\label{eq:AWW-coupling}
    \frac{g^2_2}{32\pi^2} \frac{A}{F_A} W_{\mu \nu}^i\widetilde{W}^{i\mu \nu}\,
\end{equation}
Integrating out the $SU(2)_L$ instantons, we obtain the following EW axion potential where $\bm W$ and $\bm{\widetilde{W}}$ are the electroweak gauge field and its dual~\cite{Nomura:2000yk}
\begin{equation}
   V=\frac{1}{2}\Lambda_A^4 (1-\cos(A/F_A)),
\end{equation}
where
\begin{equation}\label{eq:potential-EW-instanton-SUSY}
    \Lambda_A^4 = 2c\times \epsilon^{10} m_{3/2}^3 M_{\rm PL} e^{-\frac{2\pi}{\alpha_2(M_{\rm PL})}}.
\end{equation}
Here, $c$ is an $\mathcal{O}(1)$ constant, $\epsilon\simeq 1/17$ \cite{Buchmuller:1998zf} the flavor suppression factor due to the Froggatt-Nielsen symmetry \cite{Froggatt:1978nt} and $\alpha_2(M_{\rm PL}) \simeq 1/23$ the $SU(2)_L$ gauge coupling constant  at the Planck scale $M_{\rm PL}\simeq 2.4\times 10^{18}$ GeV. We reproduce the observed dark energy at around the hill top of the EW axion potential, that is \cite{Planck:2018vyg, Workman:2022ynf}
\begin{equation}\label{eq:CC-observed}
    \Lambda_{\rm obs}^4 \approx 7 \times 10^{-121} M_{\rm PL}^4,
\end{equation}
for $c\simeq 1$ and the gravitino mass $m_{3/2}\simeq1.8$ TeV.

However, in the non-SUSY SM, the $SU(2)_L$ gauge coupling constant is smaller as $\alpha_2(M_{\rm PL}) \simeq 1/48$, and hence the potential energy of the EW axion becomes too small. This can be partially remedied by removing the flavor suppression factor $\epsilon^{10}$. For instance, a model with a discrete $Z_{10}$ Froggatt-Nielsen symmetry can generate the axion potential without the flavor suppression factor~\cite{Choi:2019jck}.
This $Z_{10}$ Froggatt-Nielsen symmetry is anomaly free and may originate from a gauge  symmetry~\cite{Choi:2019jck}. %
Alternatively, we can simply consider cases without the flavor symmetry in the non SUSY standard model.
Then, the maximal potential energy is given by~\cite{Ibe:2018ffn}
\begin{equation}\label{eq:potential-nonSUSY-Ibe}
    \Lambda_A^4 = c'\times M_{\rm PL}^4 e^{-\frac{2\pi}{\alpha_2(M_{\rm PL})}},
\end{equation}
which still gives the \textbf{ten} orders of magnitude smaller than the observed dark energy in Eq.~\eqref{eq:CC-observed}.\footnote{It is pointed out that we can explain the observed vacuum energy if we adopt a larger cut-off scale than the Planck scale $M_{\rm PL}$ and c' a larger than $\mathcal{O}(1)$ \cite{McLerran:2012mm}.}

Now we return to the case of SM + one complex $SU(2)_L$ triplet Higgs.
From Eq.~\eqref{eq:a2beta}, it is clear that $\alpha_2^{-1}(M_{\rm PL})$ becomes larger than that in SM. The shift is estimated as $-1/(3\pi) \ln (M_{\rm PL}/1{\rm TeV}) \approx -3.8$, which gives $\exp(2\pi  (3.8)) \approx 2.3 \times 10^{10}$. Therefore, the potential energy of the EW axion explains the observation Eq.~(\ref{eq:CC-observed}) at $A\simeq \pi F_A$.

\section{Discussion}

In all of the above discussions on the potential energy of the EW axion, we have assumed that the vacuum energy of the true potential minimum vanishes. This is at least consistent with the de-Sitter conjecture on the quantum breaking of the classical gravity \cite{Dvali:2014gua} or on the string swampland \cite{Vafa:2005ui}, or our EW axion dark energy hypothesis is consistent with the anthropic principle \cite{Weinberg:2000yb}.\footnote{The dark energy at the present universe may consist of two components unless the axion sits at its potential minimum already. One is truly the cosmological constant and the other the potential energy of the EW axion. The total dark energy is then constrained by the anthropic principle \cite{Weinberg:2000yb}. If both terms are comparable and the cosmological constant is positive, we do not have a severe fine-tuning problem discussed below even if the decay constant is $F_A\simeq 10^{16}$ GeV for a good quality of the EW axion. And the cosmic birefringence can be natuarlly explained by the movement of the EW axion. However, this positive cosmological constant is inconsistent with the de-Sitter conjecture. For the model to be consistent with the de-Sitter conjecture we must consider the cosmological constant to be the vacuum energy of some other light field $B$.} But if we want to understand the vanishing cosmological constant dynamically,
we need some deeper insight of quantum gravity which is beyond the scope of this paper.\footnote{If there is some symmetry acting on the wave function of the universe, like the chiral symmetry acting on massless fermionic fields in quantum field theory, it might be possible to impose the cosmological constant to vanish.}

However, even if the quintessence axion hypothesis is interesting, we have two serious problems. One is the quality problem of the shift symmetry as quantum gravity is believed to explicitly break any global symmetries.
The other is the fine tuning problem of the initial axion field value. Namely, we have to assume that the initial value of the EW axion $A_{\rm i}$ is taken near the hill top of the potential, that is, $A_{\rm i}\simeq \pi F_A$. The initial condition would not be a serious problem if $F_A \simeq M_{\rm PL}$. However, if it is the case, the axion quality problem becomes very serious. We describe such a conflict in a more quantitative way as follows.

As for the breaking of the shift symmetry, we have probably many processes in quantum gravity. Let us discuss wormhole instanton effects here, since it can be estimated in a relatively straightforward way by solving wormhole-axion solutions \cite{Giddings:1987cg}. The detailed calculation is given in \cite{Alonso:2017avz} in which the relevant Euclidean action of the wormhole is given by
\begin{equation}\label{eq:Swh}
    S_{wh} \simeq 0.14\times \sqrt{8\pi} \frac{M_{\rm PL}}{F_A}\,.
\end{equation}
Imposing a sufficiently small suppression, i.e., $\exp(-S_{wh})<10^{-120}$, we obtain a constraint on $F_A < 0.7\times 10^{16}$ GeV.

However, we consider that the above result is speculative due to the absence of the detailed knowledge of a quantum gravity. Maybe we have some non-perturbative contributions $\gamma$ \cite{Giddings:1987cg} to the wormhole action $S_{wh}$ which may make the total wormhole action larger as long as no cancellations in different contributions take place. Thus, we consider the decay constant $F_A \simeq 10^{16}$ GeV is sufficiently small to suppress the explicit breaking of the axion shift symmetry by the wormholes.\footnote{For a recent discussion on the quality problem of axions in a linear $\sigma$ model see \cite{Choi:2022fha}.}

On the other hand, in order to keep the equation of state of the axion to be within the observational bound, i.e., $|w+1|<0.05$ \cite{Planck:2018vyg}, the initial $A_I$ needs to be extremely close to the hill top of the potential. From a WKB approximation to the dynamics of the EW axion (and combining Eq. (10) in \cite{Choi:2021aze}), one can show that
\begin{equation}\label{eq:Initial-Condition}
    \left|\frac{A_I}{F_A} -\pi \right| \le \frac{M_{\rm PL}}{F_A} e^{-\mathcal{O
    }(M_{\rm PL}/F_A)} \,.
\end{equation}
The above is $\sim10^{-4}$ for $F_A\simeq2\times10^{17}$\,GeV but is $\sim 10^{-100}$ for $F_A\simeq10^{16}$\,GeV. Thus, an $F_A$ that is sufficiently small to keep a good quality of the EW axion will need an extreme fine-tuning of the initial condition if the EW axion is the dark energy.

To avoid the above severe conflict between the EW axion quality and the naturalness of the initial condition, one can either (1) invoke some mechanism to suppress the wormhole effects for the shift symmetry breaking; (2) study some ways to alleviate or eliminate the fine tuning of the initial condition; or (3) assume the dark energy has some other origin and allow the EW axion to be a fraction of dark matter.

A solution achieving the above (1) may be given in higher dimensional space-time where the extra dimensions are compactified. In this case, the above wormhole solution in the four dimensional Euclidean space does not exist above the compactification scale $\Lambda_{com}$. We should take into account a momentum cut-off at the compactification scale  $1/\Lambda_{com}$ in the estimation of Eq. \eqref{eq:Swh}, which would lead to a large wormhole Euclidean action \cite{Giddings:1987cg}
\begin{equation}\label{eq:Swh-extradimension}
    S_{wh} \simeq \frac{3\pi^2}{8}\frac{M_{\rm PL}^2}{\Lambda_{com}^2}.
\end{equation}
Provided that $F_A\simeq \Lambda_{com}$ we obtain $F_A \leq 2.6\times 10^{17}$ GeV to have a good quality for the quintessence EW axion. Although we can not estimate the shift-symmetry breaking effect above the compactification scale, we consider $F_A \simeq 2.6\times 10^{17}$ GeV is sufficiently small to protect the quintessence axion quality. In this case, the degree of the fine-tuning of the initial condition becomes much milder, i.e., $|{A_I}/{F_A} -\pi| \le 10^{-3}$. Additionally, with $F_A\simeq 2.6\times 10^{17}$ GeV we can explain the cosmic birefringence if the EW axion is rolling now \cite{Choi:2021aze}.

Another viable way is to consider two axions, both of which couple to the strong $SU(3)_C$ and weak $SU(2)_L$ gauge fields~\cite{Kim:2004rp}. We consider the following Lagrangian:
\begin{equation}\label{eq:two-axions}
\begin{split}
\mathcal{L} \supset& \frac{g_2^2}{32\pi^2}
\left(m_1 \frac{A_1}{F_{1}}
+m_2 \frac{A_2}{F_{2}}\right) W_{\mu \nu}^i \widetilde{W}^{i\mu \nu} \\
&~~~~~+\frac{g_3^2}{32\pi^2} \left(n_1 \frac{A_1}{F_{1}}
+n_2 \frac{A_2}{F_{2}}
\right) G_{\mu \nu}^a \widetilde{G}^{a\mu \nu},
\end{split}
\end{equation}
where $m_{1,2}$ and $n_{1,2}$ are anomaly coefficients, and $\bm{G}$ is a gluon field. It is assumed that $F_1 \sim F_2 \sim 10^{16}\,{\rm GeV}$ to allow good qualities for both axions. Assuming $m_1n_2-n_1m_2\neq0$, in the following basis,
\begin{equation}
\left(
\begin{array}{c}
    A  \\
    A^\prime
\end{array}
\right)
=
\left(
\begin{array}{cc}
    \cos\beta & -\sin\beta  \\
     \sin\beta & \cos\beta
\end{array}
\right)
\left(
\begin{array}{c}
    A_1  \\
    A_2
\end{array}
\right)
\end{equation}
where $\tan\beta=\frac{n_1F_2}{n_2F_1}$, only $A'$ couples to $\bm{G} \widetilde{\bm{G}}$ and obtains a large mass through the QCD non-perturbative effects. ({$A'$ is nothing but the QCD axion, solving the strong CP problem.}) Due to the large mass, $A^\prime$ is quickly stabilized to the minimum and can be ignored. The other orthogonal component $A$ remains massless with the following effective potential\footnote{We omit the sign.}
\begin{equation}\label{eq:A-effective}
    \mathcal{L}_{\rm eff} =  \frac{g_2^2}{32\pi^2}\frac{A}{F_{\rm eff}}W_{\mu \nu}^i \widetilde{W}^{i\mu \nu}\,,
\end{equation}
with
\begin{equation}\label{eq:F_eff}
    F_{\rm eff} = \frac{\sqrt{n_1^2F_2^2+n_2^2F_1^2}}{|m_1n_2-m_2n_1|}\,.
\end{equation}
Then, $A$ only couples to $\bm{W}\widetilde{\bm{W}}$ and can be identified as our EW axion.
A larger effective decay constant $F_A = F_{\rm eff} \sim2\times10^{17}$ GeV can be obtained with $n_{1,2}\sim10$ but a small $m_1 n_2 -m_2 n_1$ (e.g., $m_1=m_2=1$, $n_2=10$ and $n_1=9$) for instance.

For possibility (3), the EW axion can starts from a general point between $0$ and $\pi F_A$. Then, the oscillation of the EW axion after $H\sim m_A$ makes
it become a small fraction of dark matter. In fact, it is encouraging that we can naturally explain the observed cosmic birefringence \cite{Minami:2020odp,Diego-Palazuelos:2022dsq,Eskilt:2022cff} by the axion dark matter oscillation \cite{Lin:2022niw}.\footnote{This is also discussed in \cite{Fujita:2020ecn,Fujita:2020aqt}. However, they assume $F_A = M_{\rm PL}$ and the obtained parameter space is different from that in \cite{Lin:2022niw}. Furthermore, they may have a serious axion quality problem with such a large $F_A$.}

\begin{acknowledgments}
 T.\ T.\ Y.\ deeply thanks Li Fu for her helpful suggestions and encouragement and he is supported in part by the China Grant for Talent Scientific Start-Up Project and by Natural Science Foundation of China (NSFC) under grant No.\ 12175134 as well as by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. N. Y. is supported by a start-up grant from Zhejiang University.
\end{acknowledgments}

\providecommand{\href}[2]{#2}\begingroup\raggedright