\section{Introduction}
Clustering is the task of grouping objects to structures (also called clusters) that maintain internal homogeneity and external separation. Objects in the same cluster should be similar to each other while objects in different clusters should not.  In past systems, the incoming data was stored on the memory and the clustering task was performed offline. However, modern systems require clustering of streaming data as well as clustering of large datasets which may be treated as streaming data.

Clustering streaming data is a very challenging problem. Unlike traditional clustering algorithms where the dataset can be stored and scanned multiple times, clustering streaming data has to satisfy constraints such as limited memory, real-time response, unknown data statistics and a dynamic number of clusters. For example, imagine the task of segmentation by clustering of a 4K image (typical resolution of $4096 \times 2160$).  Creating an affinity matrix for traditional clustering algorithms like Normalized Cuts \cite{Shi:2000:NCI:351581.351611} and Affinity Propagation  \cite{Frey07clusteringby} is impossible because  the size of the matrix, $(4096 \times  2160) \times (4096 \times  2160)$,  cannot fit the memory. Applying other traditional algorithms like  K means  \cite{Lloyd82leastsquares}, DBSCAN \cite{ester1996density} and Meanshift \cite{Fukunaga:2006:EGD:2263309.2268796} is also limited since parameters selections for those algorithms require multiple scans of the dataset - a computationally expensive task.

A good clustering algorithm for streaming data should be able to identify the underlying structure of the data in a dynamic environment. It should respond quickly to the changing statistics of the data, be robust to outliers and process the data in real-time. It should also perform as good as if it had the global view of the data \cite{toshniwal2013clustering}.

There are two approaches to process the incoming data: the incremental learning approach and the two-phase learning approach \cite{nguyen2015survey}. Assuming the data arrives in batches that is, $n_1$ data points arrive at interval $t_1$, $n_2$ at $t_2$ and so on. In the incremental learning approach, the clustering model is incrementally updated after each batch is processed. In this approach, the clustering result is available immediately after the arrival of the last batch but a batch processing is computationally expensive. Unlike the incremental learning approach, the two-phase approach divides the learning process into two phases: online and offline. In the online phase, the arriving data is summarized and only prototypes of each batch are saved to the memory. Then, in the offline phase, after the last batch arrives, clustering is performed on the prototypes of all the batches to create the final clustering model. This approach is able to cluster high-speed data streams.

In this paper, we propose a novel online clustering algorithm based on the two-phase framework. This algorithm, which is based on a Penalized Weighted Gaussian Mixtures Model (PWG), clusters the data as it streams without knowing the number of clusters apriori. By adding a penalty function on the mixing coefficients to the traditional Gaussian Mixtures Model the algorithm can estimate the number of clusters in the dataset dynamically. Moreover, by using a Gaussian Mixtures Model, the algorithm can model any probability density functions to any required level of accuracy.

One possible application of the algorithm is online segmentation of images in modern cameras. Usually, modern cameras are equipped with an Image Signal Processor (ISP) that processes the raw data from the CCD sensor. Common blocks in such ISP are demosaicing, noise reduction and image sharpening. Every block processes the raw data in parts.  Based on the proposed method,  a segmentation block with hardware and software phases can be designed.  In the first phase, the hardware segments only part of the raw input data by using the PWG model saving weighted prototypes to the DDR. In the second phase, the software clusters the weighted prototypes to create the final segmentation result.

The rest of this paper is organized as follows.  Section \ref{sec:Related Work}  surveys various online clustering methods that are based on the two-phase framework. Section \ref{sec:opwg} presents the proposed algorithm -  Online Clustering by Penalized Weighted Gaussian Mixtures Model (OPWG). Section \ref{sec:results} presents the algorithm evaluation on both synthetic and real datasets and discusses its result and section \ref{sec:conc} concludes the paper.

\section{Related Work}
\label{sec:Related Work}

Many algorithms use the two-phase framework to extend traditional clustering algorithms to streaming data. For example, the STREAM \cite{o2002streaming} algorithm is based on the classical k-median algorithm \cite{Jain:1988:ACD:46712}. It consists of two phases. Given the streaming data is divided into buckets, in the first phase, the algorithm finds $K$ clusters in each bucket by the k-median clustering algorithm. For each bucket, it stores the clusters' centroids and their corresponding weights. In the second phase, the weighted centroids are clustered into a small number of clusters, resulting in the final data partition.

Similar to the STREAM algorithm, the Online Fuzzy C-means (OFCM) \cite{hore2008online} is based on the classical algorithm Fuzzy C-means \cite{Bezdek:1981:PRF:539444}  and it also has two phases. Assuming the data arrives in chunks, in the online phase, the arriving data is being clustered by the FCM algorithm. The clustering result is a set of centroids representing each chunk and their corresponding weights. In the post-processing phase, a weighted FCM algorithm is applied to cluster the centroids considering their weights.

Although those algorithms perform well, they have some limitations: STREAM  and OFCM need to be provided with the number of clusters, $K$, which may vary from batch to batch. Choosing a wrong $K$ may affect the algorithm as following - when all data in a given chunk comes from only one class - an over-segmentation is created. Such over segmentation does not result in information loss but increases the algorithm sensitivity to noise and its convergence time. When all the data in a given chunk comes from multiple classes, an under segmentation may occur depending on the value of $K$. Another limitation of those algorithms is the creation of spherical clusters due to the objective function of those algorithms. Creating spherical clusters can cause wrong clustering result since the clusters may have arbitrary shapes.

A hierarchical approach for clustering streaming data is proposed by the CluStream algorithm \cite{Aggarwal:2003:FCE:1315451.1315460}. The CluStream algorithm extends the notion of a feature vector of the BIRCH \cite{Zhang:1996:BED:235968.233324} method to create micro-clusters in the online phase. In the offline phase, the micro-clusters are clustered into bigger clusters using the K-means algorithm. Like the above methods, the number of clusters should be provided to this algorithm.

A density-based online clustering method is the DenStream \cite{cao2006density}.  The DenStream algorithm extends the notion of core points of the  DBSCAN algorithm to a new concept of micro-clusters. It also has two phases. In the online phase, the algorithm maintains micro clusters structures which approximately capture the density of the data stream. In order to create the final data partition, in the offline phase, a variant of the DBSCAN algorithm is applied to those clusters. The DenStream algorithm should be provided with a cluster radius threshold and data fading rate. Providing the cluster radius is a main drawback of the algorithm. As stated, in streaming data, the data statistics may vary along the time and therefore a single cluster radius may not fit all the data.

\section{Online Clustering by Penalized Weighted GMM}
\label{sec:opwg}

Gaussian Mixtures Model (GMM) is usually the model of choice in distribution-based clustering.  In this approach, the data objects are assumed to be generated according to a mixture of Gaussians and the clusters are defined as objects belonging to the most likely Gaussian. Therefore, finding the clusters of a given data is equivalent to estimating the parameters of the GMM.

The classical GMM regards the importance of each object in the dataset equally and the clustering result is equally affected by all objects. More generally, different objects may have different impacts on the dataset. An object may have a weight that represents its "mass". The larger the weight the larger the dominance of the object. In such a case, the clustering result will be biased by those objects with high dominance. Such an approach takes place in the second phase of STREAM and OFCM, where the algorithms solve a weighted clustering problem.

The OPWG algorithm works similarly. Assume a streaming data arrives and is processed in batches so that $n_1$ data points arrive at time $t_1$, $n_2$ at $t_2$ and so on, OPWG clusters the data in two phases. In the first phase, the algorithm clusters each batch based on a degenerated PWG model where all the weights equal to 1, unless there is a prior knowledge on the data. Each batch is clustered to $K < K_{max}$ Gaussians mixture, whose centroids and mixing coefficients are saved. The number of clusters is found dynamically starting from $K_{max}$ and converging to a smaller value $K$. The mixing coefficient acts as a weight - the larger the coefficient - the more dominant its corresponding centroid is.
In the second phase, the algorithm clusters the weighted centroids from all the batches to get the final clustering results. As in the first phase, this phase also estimates the number of clusters. The PWG model is discussed next.

\subsection{Penalized Weighted GMM}

The OPWG algorithm is based on a penalized weighted GMM, which will be defined next. Let ${x_i} \in \mathbb{R}^{d}$ be a random vector following a multivariate Gaussian distribution with  mean $\mu \in \mathbb{R}^{d} $ and a covariance matrix $\Sigma \in \mathbb{R}^{d \times d} $. Let $w_i>0$ be a weight indicating the dominanace of an observation ${x_i}$. As suggested in \cite{gebru2015algorithms} the weights $w_i$ are incorporated into the regular GMM model by "observing ${x_i}$, $w_i$ times", resulting with the following PDF of weighted GMM:
\begin{equation}
p(x_i;\Theta ,w_i) = \sum_{k=1}^K {\pi_k N(x_i;\mu_k,\frac{1}{w_i}\Sigma_k) },
\label{eq:gmm}
\end{equation}

where $\Theta={\{\pi_1,...,\pi_K;\theta_1,...,\theta_K\}}$ are the mixtures parameters: $\{\pi_1,...,\pi_k\}$ are the mixing coefficients satisfying $\sum_1^K \pi_k =1$ and $\theta_k=\{\mu_k,\Sigma_k\}$ are the parameters of the $k^{th}$ component. $K$ is the number of components in the model and it corresponds to the number of clusters in the data. Let $X=\{x_1,...,x_N\}$ and $W=\{w_1,...,w_N\}$  be the i.i.d observed data and its corresponding weights, respectively,  the observed data log-likelihood is:
\begin{equation}
l(\Theta) = \sum_{i=1}^N {\ln(\sum_{k=1}^K {\pi_k N(x_i;\mu_k,\frac{1}{w_i}\Sigma_k)})}.
\label{eq:ll}
\end{equation}

Traditional clustering algorithms require the user to provide the number of clusters. In cases where this information is unavailable,  it is common to run the algorithm several times with different values of $K$ and to choose the one that minimizes some information criteria such as AIC \cite{akaike1998information} or BIC \cite{schwarz1978estimating}. In online clustering, such approach of model selection is not applicable.
We therefore consider a different approach, in which we start with a large number $K_{max}$ of components and shrink insignificant mixing coefficients to zero to preserve a suitable number $K$ of components. This elimination can be obtained by adding a penalty to the objective function on the mixing coefficients. It is shown in \cite{PenHuaZha16} that good results can be achieved  by choosing a penalty of the form of $\ln{\frac{\epsilon + \pi_k}{\epsilon}} = \ln{(\epsilon + \pi_k)} - \ln{\epsilon} $, where $\epsilon$ is a small positive number. This penalty is a monotonically increasing function of $\pi_k$, and it goes to zero when a mixing coefficient, $\pi_k$, goes to zero.

Considering this penalty, the penalized log-likelihood function becomes:

\begin{equation}
l_p(\Theta) = l(\Theta) - N \lambda D_f \lbrack \sum_{k=1}^K {\ln{(\epsilon + \pi_k)} - \ln{\epsilon}} \rbrack,
\label{eq:llp}
\end{equation}

where $l(\theta)$ is the log-likelihood, $\lambda$ is a tuning parameter and $D_f$ is the number of free parameters of each mixture component. Specifically, for GMM with a full covariance matrix, $D_f = \frac{d(d+1)}{2}+d+1$, and for GMM with a diagonal covariance matrix, $D_f =d+d+1$. Integrating \eqref{eq:ll} and \eqref{eq:llp} results in a penalized weighted GMM log likelihood function:

\begin{equation}
\begin{split}
l_{pw}(\Theta) &=  \sum_{i=1}^N {\ln(\sum_{k=1}^K {\pi_k N(x_i;\mu_k,\frac{1}{w_i}\Sigma_k)})}\\
&- N \lambda D_f \lbrack \sum_{k=1}^K {\ln{(\epsilon + \pi_k)} - \ln{\epsilon}} \rbrack.
\end{split}
\label{eq:llpw}
\end{equation}

In order to estimate the model parameters, $\Theta$, we maximize the penalized log-likelihood function.  However, as in mixtures model, a direct maximization of the penalized log-likelihood function is analytically intractable. We therefore use the iterative EM  algorithm \cite{dempster1977maximum} to find a solution.

\subsection{EM algorithm}

Let us denote by  $Z=\{z_1,...,z_N\}$ a set of hidden variables associated with the observed data $X=\{x_1,...,x_N\}$ such that $z_{ik}=1$ iff $x_i$ is generated by the $k^{th}$ component of the mixture and $z_{ik}=0$ otherwise. Calculating the conditional expectation of the penalized log-likelihood of the complete data given the hidden variables results in:

\begin{equation}
\begin{aligned}
&Q(\Theta,\Theta^t) = E[\ln p(X,Z|\Theta)|X,\Theta^t] = \\
&\sum_{i=1}^N {\sum_{k=1}^K {\eta_{ik}^t \big \lbrack \ln{\pi_k^t} -\ln{|\Sigma_k^t|^\frac{1}{2}}   -\frac{w_i}{2} (x_i - \mu_k^t)^T (\Sigma_k^t)^{-1}(x_i-\mu_k^t) \big \rbrack }} \\
&  -N\lambda D_f \lbrack \sum_{k=1}^K {\ln{(\epsilon + \pi_k)} - \ln{\epsilon}}  \rbrack,
\end{aligned}
\label{eq:ellpw}
\end{equation}
where the posterior probabilities are updated with:
\begin{equation}
\begin{split}
\eta_{ik}^{t+1} &= p(z_i=k|x_i;w_i,\Theta^t) \\
&= \frac{\pi_k^t N(x_i;\theta^t,w_i)}{\sum_{k=1}^K \pi_k^t  N(x_i;\theta^t,w_i)} \\
&= \frac{\pi_k^t N(x_i;\mu_k^t,\frac{1}{w_i}\Sigma_k^t)}{\sum_{k=1}^K \pi_k^t N(x_i;\mu_k^t,\frac{1}{w_i}\Sigma_k^t)}.
\end{split}
\label{eq:estep}
\end{equation}

In the M-step, the parameters $\Theta={\{\pi_1,...,\pi_K;\mu_1,...,\mu_K;\Sigma_1,...,\Sigma_K\}}$ are updated by maximizing  \eqref{eq:ellpw}:

\begin{equation}
\begin{split}
\Theta^{t+1} = arg \displaystyle\max_\Theta Q(\Theta,\Theta^t).
\end{split}
\end{equation}

The parameters $\{\pi_1,...,\pi_K\}$ and $\{\mu_1,...,\mu_K,\Sigma_1,...,\Sigma_K\}$ can be updated separately since they are not intervened at \eqref{eq:ellpw}. This results in:

\begin{equation}
\mu_k^{t+1} = \frac{\sum_{i=1}^N {w_i \eta_{ik}^{(t+1)}x_i}}{\sum_{i=1}^N {w_i \eta_{ik}^{(t+1)}}},
\label{eq:mstep_mu}
\end{equation}

\begin{equation}
\Sigma^{t+1}_k = \frac{\sum_{i=1}^N {w_i \eta_{ik}^{(t+1)}(x_i-\mu_k^{t+1})(x_i-\mu_k^{t+1})^T}}{\sum_{i=1}^N {\eta_{ik}^{(t+1)}}}.
\label{eq:mstep_sigma}
\end{equation}

Incorporating the constraints on $\pi_k$ and using Lagrange multipliers, we obtain:

\begin{equation}
\pi_k^{t+1} = \max \{0,\frac{1}{1-K\lambda D_f} [\frac{1}{N} \sum_{i=1}^N{\eta_{ik}^{t+1}} -\lambda D_f]\}.
\label{eq:mstep_pik}
\end{equation}

In summary, given predefined weights, $W$, the EM algorithm starts with $K_{max}$ which is an upper bound on the number of clusters in the data. In each iteration, a mixing coefficient may be shrunk to zero by
\eqref{eq:mstep_pik}, eliminating its corresponding component.  In this way, in each iteration, the algorithm starts with a number of components and may converge to a smaller number. Note that when all the weights are queal to 1's, the algorithm converges to the algorithm proposed by \cite{PenHuaZha16} and will be referred to as Penalized Gaussian Mixtures Model (PGMM).

\subsubsection{Parameters selection}

The PWG algorithm requires three parameters: the tuning parameter $\lambda$, the initial number of clusters $K_{max}$, and the covariance matrix type.

Equation \eqref{eq:mstep_pik} constrains the choice of the tuning parameter and the number of clusters. In this equation, we want a cluster to be eliminated when its average responsibilities, $\frac{1}{N} \sum_{i=1}^N{\eta_{ik}^{t+1}}$,  is less than the threshold $\lambda D_f$. Therefore,  the term $\frac{1}{1-K_{max}\lambda D_f}$ should be positive. In addition, $\lambda$ should be non-negative from its definition. Therefore :

\begin{equation}
  0  \leq \lambda < \frac{1}{K_{max}D_f}.
\end{equation}

 When $\lambda$ is close to its upper limit, the algorithm will prefer a small number of components. When $\lambda$ approaches zero, the algorithm will converge to the regular Weighted GMM (WGMM) algorithm, proposed by \cite{gebru2015algorithms}. In \cite{PenHuaZha16} it is suggested to choose $\lambda$ by running the algorithm with different values and choosing the one that minimizes the BIC criterion \cite{schwarz1978estimating}:

\begin{equation}
\lambda = \min_{\lambda} BIC(\lambda).
\label{eq:lambda_BIC}
\end{equation}

While the choice of the tuning parameter is constrained, the choice of the type of the covariance matrices is not. Thus, the covariance matrices can be either full or diagonal. However, the best practice is to choose diagonal covariance matrices due to three reasons. First, modeling of a full covariance GMM can be equally achieved by using diagonal covariance GMM with more components. Second, a GMM with diagonal covariance matrices is more computationally efficient. Third, the results of a GMM with diagonal matrices are empirically superior to GMM with full covariance matrices \cite{Bimbot:2004:TTS:1289340.1289376}.

\section{Algorithm Evaluation}
\label{sec:results}

The proposed method has been evaluated by several experiments on both synthetic and real datasets. The synthetic datasets included: a) the data used by scikit package \cite{scikit-learn}  in its clustering algorithm comparison b) simulated Gaussian mixtures PDFs with a different number of components.
For real datasets evaluation, we used OPWG to segments images from Berkeley BSD dataset \cite{MartinFTM01}.

\subsection{Algorithm Evaluation - Synthetic Datasets}

To evaluate the algorithm on synthetic data we performed 50 experiments on each dataset and averaged the results of the F1 score and the Normalized Mutual Information (NMI) metric, as defined as in \cite{Pfitzner:2009:CES:1554488.1554491}.

Every dataset consists of $N=100,000$ 2D data samples which is further divided into equal batches of $n_i=1000$ samples each.
Since the order matters when clustering streaming data, two modes of operations were considered: Mode A in which the data arrives randomly and each batch may contain samples from all the clusters in the dataset and Mode B in which the data in each batch arrives from a subset of the dataset and may not contain samples from all the underlying clusters. To simulate the latter case, the 2D data has been sorted according to x (or y) coordinate and then divided into batches.

The tuning parameter $\lambda$ was chosen to be $0.005$ in the online stage. Different values of $\lambda$ in the range $[0.003-0.006]$ were tested in the post-processing stage and the one that minimized the BIC of the model was chosen. The covariance matrices were chosen to be diagonal in the online stage and full in the post-processing stage. The maximum number of clusters was chosen to be $K_{max}=25$.

In our evaluation, we compared our proposed algorithm to the following algorithms: OFCM  with a number of clusters equal to that found by OPWG and a batch size of $n_i=1000$ samples, PGMM running on the whole dataset with $\lambda=0.005$, $K_{max}=25$ and a full covariance matrix, and GMM  running on the whole dataset with a number of clusters provided by PGMM and with a full covariance matrix.

The results for mode A are reported in Table \ref{table:6db_random} and illustrated in Figure \ref{fig:6db_random}.
As expected, when the full stream is available, as in PGMM and GMM, better results are achieved, whereas, in online mode, the performance deteriorates.
The underlying assumption in online mode that all batches are independent and have an identical distribution yields sub-optimal solution. Further, comparing OPWG with OFCM shows that except for DB2, OPWG outperforms OFCM although its initial number of clusters $K_{max}$ is bigger. In this mode, the average number of clusters that OPWG found, $K=3.18, 1.92, 2.96, 2.48, 2,1 $, was close to the number of clusters that PWG  found - $K= 3.34,5.22,3.04,2.38,2.14,2.76$.

\begin{table*}[h]
\begin{center}
\begin{tabular}{lllllllllll}
\hline
            & \multicolumn{2}{c}{DB1} & \multicolumn{2}{c}{DB2}       & \multicolumn{2}{c}{DB3}       & \multicolumn{2}{c}{DB4}       & \multicolumn{2}{c}{DB5}       \\
            & F1              & NMI   & F1            & NMI           & F1            & NMI           & F1            & NMI           & F1            & NMI           \\ \cline{2-11}
OPWG  & \textbf{0.46}   & 0.00  & 0.54          & 0.00          & \textbf{0.98} & \textbf{0.92} & 0.75          & 0.63          & 0.77          & 0.73          \\
OFCM & 0.40            & 0.00  & \textbf{0.75} & 0.19          & 0.91          & 0.78          & 0.78          & 0.60          & 0.77          & 0.73          \\
PGMM         & 0.38            & 0.00  & 0.57          & \textbf{0.55} & \textbf{0.98} & \textbf{0.92} & 0.86          & \textbf{0.82} & 0.79          & 0.74          \\
GMM         & 0.38            & 0.00     & 0.56          & 0.54          & \textbf{0.98} & \textbf{0.92} & \textbf{0.85} & 0.81          & \textbf{0.80} & \textbf{0.77} \\ \hline
\end{tabular}
\caption
{Mode A. Clustering quality measures for each of the algorithms discussed for all datasets considered apart from the last one, which is noise.}\label{table:6db_random}
\end{center}
\end{table*}

\begin{figure}[htb]
\centering
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__1_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__1_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__1_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__1_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
\begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__2_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__2_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__2_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__2_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__3_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__3_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__3_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__3_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__4_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__4_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__4_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__4_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__5_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__5_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__5_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__5_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__6_opwg-eps-converted-to.pdf}
    \caption{OPWG}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__6_ofcm-eps-converted-to.pdf}
    \caption{OFCM}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__6_pgmm-eps-converted-to.pdf}
    \caption{PGMM}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/random_1__6_gmm-eps-converted-to.pdf}
    \caption{GMM}
  \end{subfigure}\\
\caption
{Mode A.  Each algorithm was evaluated on six different datasets. For OPWG and PGMM, the initial number of clusters was 25. OFCM and GMM were initialized with the number of clusters found in OPWG and PGMM, respectively.}\label{fig:6db_random}
\end{figure}

The results for mode B are reported in Table \ref{table:6db_sorted} and illustrated in Figure \ref{fig:6db_sorted}.  In Mode B not all clusters are represented in each batch. While this fact has a little effect on the performance of PGMM and GMM, it significantly affects the online algorithms. However, as can be seen, the OPWG is less affected than the OFCM by this phenomena. Note that, in this mode,  because only subset of the data is represented in each batch, the average number of clusters that the OPWG and PGMM found increased compared to mode A. The average number that OPWG found was $K= 4.84, 7.84, 4.48, 3.04, 2.14, 1.02 $ and the average number of clusters that PGMM found was $K=3.25, 5.4, 3.08, 2.48, 2.32, 2.9$.

\begin{table*}
\begin{center}
\begin{tabular}{lllllllllll}
\hline
            & \multicolumn{2}{c}{DB1}       & \multicolumn{2}{c}{DB2}       & \multicolumn{2}{c}{DB3}       & \multicolumn{2}{c}{DB4}       & \multicolumn{2}{c}{DB5}       \\
            & F1            & NMI           & F1            & NMI           & F1            & NMI           & F1            & NMI           & F1            & NMI           \\ \cline{2-11}
OPWG  & \textbf{0.41} & \textbf{0.10} & 0.46          & 0.51          & 0.93          & 0.83          & \textbf{0.99} & \textbf{0.98} & 0.78          & 0.73          \\
OFCM & 0.31          & 0.00          & 0.45          & 0.47          & 0.83          & 0.75          & 0.83          & 0.60          & 0.80          & 0.77          \\
PGMM         & 0.39          & 0.00          & 0.54          & \textbf{0.55} & \textbf{0.98} & \textbf{0.93} & 0.87          & 0.84          & 0.82          & 0.78          \\
GMM         & 0.38          & 0.00          & \textbf{0.55} & 0.54          & \textbf{0.98} & \textbf{0.93} & 0.86          & 0.82          & \textbf{0.83} & \textbf{0.79} \\ \hline
\end{tabular}
\caption
{Mode B. Clustering quality measures for each of the algorithms discussed for all datasets considered apart from the last one, which is noise.
}\label{table:6db_sorted}
\end{center}
\end{table*}

\begin{figure}[htb]
\centering
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_1_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_1_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_1_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_1_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_2_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_2_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_2_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_2_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
    \begin{subfigure}[b]{.17\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_3_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_3_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_3_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_3_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_4_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_4_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_4_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_4_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_5_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_5_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_5_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_5_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_6_opwg-eps-converted-to.pdf}
    \caption{OPWG}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_6_ofcm-eps-converted-to.pdf}
    \caption{OFCM}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_6_pgmm-eps-converted-to.pdf}
    \caption{PGMM}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt/sorted_6_gmm-eps-converted-to.pdf}
    \caption{GMM}
  \end{subfigure}\\
\caption
{Mode B. Each algorithm was evaluated on six different datasets. For OPWG and PGMM, the initial number of clusters was 25. OFCM and GMM were initialized with the number of clusters found in OPWG and PGMM, respectively.
}\label{fig:6db_sorted}
\end{figure}

The algorithms were further evaluated on 2D Gaussian mixtures datasets with a different number of components ($K = \{ 2, 5, 7\}$). As before, each result reported represents the average performance of running the algorithm 50 times. In this case, $\lambda$ was chosen to be $0.006$ for both modes and in mode B the dataset was arbitrarily sorted according to its y coordinate.

The results for mode A are reported in Table \ref{table:gmm_random} and illustrated in Figure \ref{fig:gmm_random}. The results for mode B are reported in Table \ref{table:gmm_sorted} and illustrated in Figure \ref{fig:gmm_sorted}. For mode A, OPWG performs well comparing to the other algorithms. Its estimation of the number of clusters, $K= 2.32, 3.92, 4.2 $, is close to the estimation of PGMM that estimated $K=2.02, 4.00, 4.44$. Moreover, although it was initialized with a bigger number of clusters, $K_{max}=25$, it performs equal or better compared to the OFCM.

Surprisingly, for mode B, the average results of  OPWG are better than for mode A. In addition, in some cases, it even performs better than the algorithms that run of the full stream. A possible reason for this is that when running in online mode, the adaptive tendency of the algorithm makes it more robust to outliers compared to running on the full stream. Thus, the algorithm achieves more accurate clustering results. For this mode, the average number of clusters that OPWG found, $K=2.48, 5.1, 6 $, is again close to the average number of clusters that PGMM found - $K=2, 4.02, 4.54$.

\begin{table}
\begin{center}
\begin{tabular}{lllllll}
\hline
            & \multicolumn{2}{c}{K=2}       & \multicolumn{2}{c}{K=5}       & \multicolumn{2}{c}{K=7}         \\
            & F1            & NMI           & F1            & NMI           & F1            & NMI              \\ \cline{2-7}
OPWG  & 0.97          & 0.90          & 0.80          & 0.76          & 0.69          & 0.69           \\
OFCM & 0.92          & 0.84          & 0.81          & 0.78          & 0.73          & 0.75           \\
PGMM         & \textbf{0.99} & \textbf{0.95} & \textbf{0.90} & \textbf{0.88} & \textbf{0.78} & \textbf{0.80}  \\
GMM         & \textbf{0.99} & \textbf{0.95} & 0.88          & 0.87          & 0.77          & 0.79          \\ \hline
\end{tabular}
\caption {Mode A - numerical results for GMM datasets. Clustering quality measures (F1 and NMI) for each of the four algorithms where each batch contains random samples from the whole dataset.}\label{table:gmm_random}
\end{center}
\end{table}

\begin{figure}[htb]
\centering
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_1_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_1_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_1_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_1_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_2_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_2_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_2_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_2_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_3_opwg-eps-converted-to.pdf}
    \caption{OPWG}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_3_ofcm-eps-converted-to.pdf}
    \caption{OFCM}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_3_pgmm-eps-converted-to.pdf}
    \caption{PGMM}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/random_3_gmm-eps-converted-to.pdf}
    \caption{GMMM}
  \end{subfigure}\\
\caption
{Mode A - results for GMM datasets. Each algorithm was evaluated on three different datasets. For OPWG and PGMM, the initial number of clusters was 25. OFCM and GMM were initialized with the number of clusters found by OPWG and PGMM, respectively.}\label{fig:gmm_random}
\end{figure}

\begin{table}[h!]
\begin{center}
\begin{tabular}{lllllllllll}
\hline
            & \multicolumn{2}{c}{K=2}       & \multicolumn{2}{c}{K=5}       & \multicolumn{2}{c}{K=7}          \\
            & F1            & NMI           & F1            & NMI           & F1            & NMI               \\ \cline{2-7}
OPWG  & 0.94          & 0.85          & \textbf{0.93} & \textbf{0.88} & \textbf{0.86} & \textbf{0.82}    \\
OFCM & 0.88          & 0.79          & 0.86          & 0.83          & 0.79          & 0.79            \\
PGMM         & 0.97          & 0.91          & 0.89          & \textbf{0.88} & 0.79          & \textbf{0.82}   \\
GMM         & \textbf{0.99} & \textbf{0.93} & 0.87          & 0.87          & 0.79          & \textbf{0.82} \\ \hline
\end{tabular}
\caption {Mode B numerical results for GMM datasets. Clustering quality measures (F1 and NMI) for each of the four algorithms where not all clusters are represented in each batch.}\label{table:gmm_sorted}
\end{center}
\end{table}

\begin{figure}[htb]
\centering
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_1_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_1_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_1_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_1_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_2_opwg-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_2_ofcm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_2_pgmm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_2_gmm-eps-converted-to.pdf}
  \end{subfigure}\\
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_3_opwg-eps-converted-to.pdf}
    \caption{OPWG}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_3_ofcm-eps-converted-to.pdf}
    \caption{OFCM}
  \end{subfigure
  \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_3_pgmm-eps-converted-to.pdf}
    \caption{PGMM}
  \end{subfigure
    \begin{subfigure}[b]{.25\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/synt_gmm/sort_3_gmm-eps-converted-to.pdf}
    \caption{GMM}
  \end{subfigure}\\
\caption
{Mode B - results for GMM datasets. Each algorithm was run on three different datasets. For OPWG and PWG, the initial number of clusters was 25. OFCM and GMM were initialized with the number of clusters found by OPWG and PGMM, respectively.}\label{fig:gmm_sorted}
\end{figure}

\subsection{Algorithm Evaluation - Real Datasets}

We tested  OPWG on real datasets for the task of image segmentation.
For this purpose, each image was converted to LAB color space and has been clustered based only on its LAB color features, i.e., segmentation by color only.
Each batch consisted of four non-overlapping rows of the image, corresponding to reading four rows from a CCD image sensor.
Minimal image size was $321 \times 482$. The algorithm was initialized with $K_{max}=25$ and the tuning parameter $\lambda$ was fixed to $0.03$ for the online stage and the best of $[0.006, 0.005, 0.004]$ for the post-processing stage. For both the online and the post-processing stages, we ran the algorithm with diagonal covariance matrices. This selection of diagonal covariance matrices is very important since full covariance matrices will result in large clusters or under segmentation. Our preference would be to do over segmentation and we thus selected diagonal matrices.

Image segmentation is not a well-defined problem and different people may segment the images differently. Thus, OPWG was not evaluated numerically as in the synthetic datasets but rather compared to the popular segmentation algorithm Efficient Graph-Based Image Segmentation (EGBIS) \cite{Felzenszwalb:2004:EGI:981793.981796}.

Figure \ref{fig:image_sementation} shows the evaluation results. As can be seen, even though the OPWG is solely based on the LAB features, it performs pretty well in segmenting the images. Comparing to the OFCM algorithm, the OPWG algorithm shows equal or better results although it was initialized with a bigger number of clusters -  $K_{max}=25$.   Even though it uses only four non-overlapping rows at each batch, the OPWG algorithm is able to achieve similar results to those of PGMM, which uses the full dataset. The OPWG algorithm sometimes even outperforms the EGBIS algorithm.

\section{Conclusions}
\label{sec:conc}
This paper proposes a penalized weighted GMM algorithm for online clustering. The algorithm is able to identify the underlying structure of the data and to estimate the number of clusters in a dynamic environment. Although the algorithm is initialized with an arbitrarily large number of clusters, results on both synthetic and real datasets show that the proposed algorithm produces partitions that are close to those you would get if you clustered the whole data at once.

Although the algorithm performs well, there is still work to be done. Better results may be achieved using different penalty functions. Moreover, since the algorithm was evaluated on datasets with low dimensionality, the effect of the "curse of dimensionality" on the algorithm is unknown. In addition, an extension of the algorithm to handle the property of data fading of streaming data should be considered.  Other directions maybe include an extension of the algorithm to deep learning architecture using RNN.

\begin{figure*}[htb]
\centering
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/bear_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/bear_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/bear_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/bear_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/bear_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/dag_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/dag_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/dag_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/dag_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/dag_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/orev_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/orev_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/orev_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/orev_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/orev_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/star_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/star_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/star_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/star_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/star_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

      \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/man_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/man_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/man_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/man_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/man_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

       \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/elephent_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/elephent_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/elephent_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/elephent_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/elephent_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

       \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/snake_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/snake_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/snake_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/snake_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/snake_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

       \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/boat_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/boat_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/boat_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/boat_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/boat_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/erez1_rgb-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/erez1_opwg-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/erez1_ofcm-eps-converted-to.pdf}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/erez1_pgmm-eps-converted-to.pdf}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/erez1_egbis-eps-converted-to.pdf}
  \end{subfigure}\\%

         \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/cat_rgb-eps-converted-to.pdf}
    \caption{Image (RGB)}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/cat_opwg-eps-converted-to.pdf}
    \caption{OPWG }\label{fig:1c}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/cat_ofcm-eps-converted-to.pdf}
    \caption{OFCM }\label{fig:1c}
  \end{subfigure
    \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/cat_pgmm-eps-converted-to.pdf}
    \caption{PGMM }\label{fig:1c}
  \end{subfigure
  \begin{subfigure}[b]{.20\linewidth}
    \centering
    \includegraphics[width=.99\textwidth]{images/seg2/cat_egbis-eps-converted-to.pdf}
    \caption{EGBIS }\label{fig:1d}
  \end{subfigure}\\%
  \caption{Color segmentation results for different images from BSD dataset. Note that the last two images are not from BSD dataset. Those images were taken by a mobile phone camera. Their resolution is $1600 \times 1200$ }
  \label{fig:image_sementation}
\end{figure*}

\pagebreak

{\small
\bibliographystyle{ieee}


\section{Introduction}
The study of various modes of excitations and the associated evolution of nuclear shapes along spin and isospin axes in atomic nuclei is one of the fundamental quests in nuclear physics. While nuclei with ``magic numbers'' of protons and/or neutrons have spherical ground states, as one moves away, the polarizing effect of added nucleons leads to deformation. Throughout the nuclear landscape, this onset of deformation is usually a gradual process, however in neutron rich nuclei around mass $A\sim100$ the shape change is rather drastic and abrupt.
The ground states of Sr and Zr isotopes with $N$ ranging from the magic number $N = 50$ up to $N<60$ are weakly deformed, however,
they undergo a rapid shape transition from nearly spherical to well deformed prolate deformations as $N = 60$ is approached.
The sudden nature of  shape transition in Sr and Zr isotopes is evident from the abrupt changes in the two neutron separation energies \cite{Hager06}
and mean-square charge radii \cite{Buchinger,Campbell}, but also from the excitation energies of 2$^+_1$ states and $B(E2)$ values \cite{nndc}.
On the other hand, in
isotopes with $Z\geq42$ the shape change is rather gradual \cite{Charlwood, Hager06} showing also characteristic signatures of triaxiality.
This strong dependence of the observed spectroscopic properties, both on the number of protons and neutrons, makes the neutron-rich $A\sim100$
nuclei an excellent mass region for testing various theoretical models.

Many
experimental and theoretical studies have already been reported on the structure of these nuclei.
More specifically for
the Zr isotopes, the onset of deformation at $N=60$ has been described by a number of theoretical models \cite{GCM,Rodriguez2010,Delaroche,Reinhard99,Skalski,sm2,sm3,sm4,Petrovici2012,Ozen2006,Garcia2005,VanIsacker,Nomuraibm,Xiang2012},
however, none of the models have been able to successfully reproduce the aforementioned rapid change.
Very recently, the abrupt shape changes were correctly described by large-scale Monte-Carlo Shell Model (MCSM) calculations
\cite{Otsuka, Togashi}.
In the so-called type-II shell evolution scenario,
the (prolate) deformed states in the isotopes with $N\geq60$ are associated with
proton excitations to
the $0g_{9/2}$ orbital.
Driven by the central and tensor components of the effective (proton-neutron) interactions,
these excitations result in a lowering and subsequent filling of the neutron $0g_{7/2}$ and $0h_{11/2}$ orbitals \cite{Togashi}.
Therefore, both protons and neutrons act coherently to induce the deformation.
These multiparticle-multihole excitation based prolate deformed  states are expected to coexist as excited 0$^+$ state along with the spherical ground state in Zr isotopes with $N<60$.
The crossing of these two distinct coexisting quantum configurations at $N = 60$ manifests as the abrupt change in structure of the ground state and has been interpreted as a quantum phase transition (QPT) from a spherical to a deformed phase \cite{Togashi}.

Appearance of low-lying $0_2^+$ states in Sr and Zr isotopes with $N<60$ supports the shape-coexistence scenario in these isotopes as an explanation of the rapid shape evolution. However, in order to have a complete picture of the phenomena it is essential to have a precise information on the deformation as inferred from the electromagnetic transition strengths for nuclei in the region of the shape transition.
Experimental programs pursued in the past to determine these parameters in the $A\sim100$ region nuclei have been able to furnish useful insight into the phenomena \cite{Clement,ClementPRC,Chakraborty2013,Kremer}.  For example, coexistence of highly deformed prolate and spherical configurations have already been established in $^{96,98}$Sr \cite{Clement,ClementPRC} and $^{94,96}$Zr \cite{Chakraborty2013,Kremer}.

Shape-coexistence has also been suggested in $^{98}$Zr \cite{Wu2004,heyde2011}, a key nucleus for understanding the QPT phenomenon in the Zr isotopes. Several attempts have been made in the past to determine the electromagnetic transition rates in $^{98}$Zr. In previous work by Bettermann {\em et al}  \cite{Bettermann}, using the $\beta\gamma\gamma$ fast-timing method,
a lifetime of $\tau$(4$^+$) = 29 $\pm$9 ps
was obtained for the 4$^+$ state in $^{98}$Zr, while only upper limits could be obtained for the 2$^+$ state  ($\tau$(2$^+_1$)$\leq$ 15 ps) and the
6$^+$ state ($\tau$(6$^+_1$) $\leq$ 14 ps).
More recently, in the work by Ansari {\em et al.} \cite{Ansari}, the
so far most precise
upper limit for the lifetime of $2^+_1$ state in  $^{98}$Zr was
obtained as $\tau(2^+_1)$ $\leq$ 6 ps, and an upper limit of 15 ps was reported for the 4$^+_1$ state.
These limits for the 2$^+$ state did not allow to prove whether the onset of collectivity at $N=60$, i.e. in $^{100}$Zr, is as rapid as expected from the drop in excitation energy,
and predicted by the MCSM, or whether the B(E2) value already increases in $^{98}$Zr. Moreover, in the absence of precise information on $B(E2)$ values one cannot draw any definite conclusion about the properties of the coexisting structures.

The present letter reports the results of lifetime measurement of the first excited 2$^+$, 4$^+$ and 6$^+$ states in $^{98}$Zr using the \RDDS\ method and confirms for the first time shape-coexistence in this nucleus. In addition, our result shows very small B($E2;2_1^+\rightarrow0_1^+$) value in $^{98}$Zr thereby confirming the very sudden onset of collectivity at $N=60$.

The excited states in $^{98}$Zr were produced in fission reactions, where a $^{238}$U beam with an energy of 6.2 MeV/u delivered by the GANIL facility,
impinged on a  2.3 mg/cm$^{2}$ thick $^{9}$Be target. The VAMOS++ \cite{Pullanhiotan2008,Rejmund2011} spectrometer positioned at 20$^\circ$ with respect to the beam axis was used to detect and identify the reaction products in mass, charge, and atomic number on an even-by-event basis.
Gamma rays emitted from the reaction products were detected using the EXOGAM array \cite{Simpson2000} consisting of ten segmented Clover detectors arranged in two rings (three detectors at 135$^\circ$ and seven detectors at 90$^\circ$) with respect to the spectrometer axis.
Correlation between the prompt $\gamma$ rays detected at the target position and the detection of an ion in the focal plane of VAMOS++ triggered the event building, thus allowing the selection of prompt $\gamma$ rays from isotopically identified fission fragments.

\begin{figure
\centering
 \includegraphics[clip=true,width=.500\textwidth, angle=0]{spectra_shifted.eps}
\caption {\label{fig2} $\gamma$-ray spectra showing the $2^+\rightarrow0^+$, $4^+\rightarrow2^+$, $6^+\rightarrow4^+$ and $3^-\rightarrow2^+$ transitions in $^{98}$Zr, observed in the detectors at 135$^\circ$ for seven target-degrader distances. The spectra are Doppler corrected using the velocity measured in VAMOS++ after the degrader. The dotted and dashed lines indicate the positions of the shifted and unshifted components, respectively. Also displayed are the gaussian fits to the spectra (green), and parameters chosen for fitting the shifted (red), and unshifted (blue) components.
}
\end{figure}
A compact plunger device was used for the lifetime measurements using the
\RDDS~(RDDS) technique
\cite{Dewald2012}. The recoiling $^{98}$Zr nuclei had an average velocity of 37.3 $\mu$m/ps and were slowed down in a 4.9 mg/cm$^2$ thick Mg degrader foil placed behind the target.  Data was collected for seven distances ranging from 37 to 1554 $\mu m$ with average running times of approximately 24 hours per distance.

The energy of the $\gamma$ rays was Doppler corrected event-by-event using the measured velocity and direction of the fission fragments as measured after the degrader foil in VAMOS++.
Figure \ref{fig2}~shows the Doppler-corrected $\gamma$-ray spectra of $^{98}$Zr observed with the EXOGAM detectors located at 135$^\circ$  for different target-degrader distances. The $\gamma$ rays that were emitted after passing through the degrader appear at the correct transition energy in the Doppler corrected spectra, whereas those emitted between target and degrader are shifted to lower energies when observed at backward angles. Decay curves $I^u_{ij}(x)$/(${I^u_{ij}(x)+I^s_{ij}}(x)$) were constructed from the intensities of the shifted ($I^s_{ij}(x)$) and unshifted ($I^u_{ij}(x)$) components of the transitions as a function of target-degrader distance. For each distance, the normalisation factors (${I^u_{ij}(x)+I^s_{ij}}(x)$) were found to be consistent with the number of ions identified in VAMOS++. This allowed us to use the relative number of ions identified in VAMOS as normalization factor for the weak and contaminated  583 keV transition depopulating the $3^-\rightarrow2^+$ states.
The lifetimes of individual states were determined from the relative intensities of the two components as a function of target-degrader distance
using the differential decay curve method (DDCM) analysis for singles plunger data \cite{Dewald1989}. All observed feeders have been taken into account as per their experimentally observed intensities. The unobserved feeding contribution was subtracted following the procedure explained in Ref. \cite{Dewald1989}. We have also explored independently the errors in lifetime due to relativistic effects, possible deorientation effects  and solid angle effects, however, error due to these effects were found negligible as compared to the experimental uncertainties and are not taken into account in present analysis.
Further details on the experiment and the analysis techniques
can be found in Refs. \cite{Hagen2018,Hagen2017}.

As most important result the lifetime
of the 2$^+_1$ state was determined
for the first time
as $\tau=3.8\pm$0.8 ps.
In addition, we were also able to determine the lifetime for the 4$^+_1$ and  the 6$^+_1$ state as 7.5$\pm$1.5 ps and 2.6$\pm$0.9 ps, respectively. The $B(E2)$ values from the present work are listed in table I.
Our new lifetime measurement
yields a $B(E2;2_1^+\rightarrow0_1^+)$ value of 2.9(0.6) W.u, thereby confirming the rapid onset of deformation in Zr isotopes only at $N=60$.

\begin{table}
\caption{\label{table1}Summary of results for the lifetime measurements in $^{98}$Zr.}

\begin{ruledtabular}
\begin{tabular}{ccccccccc}

$I^\pi_i\rightarrow I^\pi_f$ &$\tau_{exp}$ (ps) &$\tau_{lit.}$ (ps)&$B(E2,\downarrow)_{exp} (W.u.)$ \\

\hline
$2^+_1\rightarrow0^+_1$     & 3.8(0.8)&$\leqslant6$\cite{Ansari}  & 2.9(0.6)$^a$(0.2)$^b$ \\
$2^+_1\rightarrow0^+_2$     &         &  & 28.3(6.0)$^a$(2.4)$^b$ \\
$4^+_1\rightarrow2^+_1$     & 7.5(1.5)&$\leqslant15$\cite{Ansari} & 43.3(8.7)$^a$(10.8)$^b$ \\
$4^+_1\rightarrow2^+_2$     &         & & 67.5(13.5)$^a$(16.9)$^b$ \\
$6^+_1\rightarrow4^+_1$     & 2.6(0.9)&   & 103.0(35.7)$^a$ \\

\end{tabular}
\end{ruledtabular}
\begin{flushleft}
\footnotesize{The errors in the $B(E2)$ values  from the measured lifetimes$^a$ and from the branching ratios$^b$ taken from \cite{Urban2017} are mentioned separately.\\
}
\end{flushleft}
\end{table}

\begin{figure
\centering
 \includegraphics[clip=true,width=.280\textwidth, angle=-90]{theory20.eps}
\caption {\label{fig4}(Color online) Comparison of calculated $B(E2;2_1^+\rightarrow0_1^+)$ values with experimental data for the Zr isotopes. Experimental data points are from refs. \cite{Raman2001,Chakraborty2013,Kumbartzki2003, Ansari,Browne2015plb}. Theoretical values are taken from refs. \cite{Nomuraibm,Mei,sm3,Togashi,Delaroche}.
}
\end{figure}

The onset of deformation in the Zr isotopes has been described by a number of theoretical models.
Information from $B(E2)$ values on the rapid structural changes in neutron-rich Sr and Zr isotopes is available from calculations
solving a five-dimensional collective Hamiltonian (5DCH) with parameters determined
within the relativistic mean-field approach
(PC-PK1 force) \cite{Mei}, the non-relativistic Hartree-Fock approach using the SLy4 \cite{Mei}
and the Gogny-D1S \cite{Delaroche} force. All these calculations suggest a picture of spherical-oblate-prolate shape transition in the neutron-rich Sr and Zr isotopes.

\begin{figure*

\centering\includegraphics[clip=true,width=0.310\textwidth,angle=-90]{98Zrbe2_ls.eps}

\caption{\label{be2_ls} Comparison between the experimental and theoretical level schemes for $^{98}$Zr. The excitation energy and spin are mentioned against each state. Labels and widths of the arrows represent calculated \cite{Delaroche,Togashi} and measured $B(E2)$ values in W.u.}

\end{figure*}

Figure \ref{fig4} compares the $B(E2;2_1^+\rightarrow0_1^+)$ values from these models with the experimental data. It can be clearly noted that the models based on the mean-field approach predict a gradual increase in collectivity already starting at $N\sim54$.
This is also the case for the interacting boson model (IBM) calculations using self-consistent mean-field approximation based on the Gogny-D1M energy density functional \cite{Nomuraibm}. However, certain differences can be noted depending on the underlying effective force. The PC-PK1 and Gogny D1S force predict a very gradual change in most of the collective properties between $N=54$ and $N=60$.
Calculations using SLy4 force, on the other hand,  predict a more pronounced change in collectivity at $N=60$ \cite{Mei}.
None of them is, however, capable to reproduce the experimental data.
This could in part be related to the fact that the particle number projection in these methods is only valid on average.
Simply speaking there might still be components in the theoretical wave functions from the neighboring isotopes with $N\pm2$, which would smoothen out
the evolution with neutron number.
It is also worth noting that the conventional large-scale shell model (SM) calculations \cite{sm3} show a rather accurate reproduction of experimental data up to $N=58$, but are constrained by valence space limitation to $N<60$.
Finally, the dramatic increase of $B(E2;2_1^+\rightarrow0_1^+)$ values between $N=58$ and 60 is very well reproduced by recent state-of-the-art MCSM calculations
\cite{Togashi}.

Using the most recent branching ratios determined by Urban {\em et al.} in the $\beta$ decay of $^{98}$Y \cite{Urban2017} we were also able to determine the $B(E2)$ values for the decay to the non-yrast states, which give important new information. The experimental excitation energies of the states and the transition strengths between them are compared with the calculations in Fig. \ref{be2_ls}.
A collective value of B($E2;2_1^+\rightarrow0_2^+$)=28.3(6.4) W.u suggests a common moderately deformed structure for the $0^+_2$ and $2^+_1$ states. Assuming a rigid, axial symmetric deformed shape yields a quadrupole deformation parameter $\beta_2\approx0.21$. These results are in general agreement with the results of MCSM calculations which predict an enhanced B($E2;2_1^+\rightarrow0_2^+$) and very small $B(E2;2_1^+\rightarrow0_1^+)$ transition strengths. However, while the calculations suggest a gradual increase in deformation of the $0^+_2$ states between $^{94}$Zr and $^{98}$Zr \cite{Otsuka}, our present results indicate a deformation which is smaller than that observed in $^{96}$Zr \cite{Kremer}.

\begin{figure

\centering\includegraphics[clip=true,width=0.280\textwidth,angle=-90]{proton_occ.eps}\\
\centering\includegraphics[clip=true,width=0.280\textwidth,angle=-90]{neutron_occ.eps}

\caption{\label{occp}(Color online) Occupation numbers of (a) protons and (b) neutrons for the states in $^{98}$Zr}

\end{figure}

As shown in Fig.~\ref{be2_ls} the $B(E2)$ values in $^{98}$Zr indicate the coexistence of three different structures at low spin:
A nearly spherical
$0^+_1$ ground state, a moderately deformed excited $0^+_2$ state as discussed above and a well-deformed (band-like) structure
possibly based on the $0^+_3$ state.
The existence of latter is supported by the rather large $B(E2)$ values within this band as shown in Fig. \ref{be2_ls}.
A significant mixing of this well-deformed structure with the moderately deformed configuration, based on the $0^+_2$ state, is indicated by the large  B($E2;4_1^+\rightarrow2_1^+$) value. A strong mixing between the two $2^+$ states was already suggested in ref.~\cite{Wu2004} based on their similar feeding and decay patterns.
Finally, the large electric monopole transition strength observed between the $0^+_3$ and $0^+_2$ states also support the mixing of two coexisting shapes \cite{heyde2011}.

The proposed triple shape coexistence scenario in $^{98}$Zr is also supported by MCSM calculations, which predict coexistence of three shapes in this nucleus. A spherical $0^+_1$ state, a prolate deformed $0^+_2$ state and a triaxial $0^+_3$ state. The underlying mechanism for stabilizing such coexisting structures has been discussed in terms of type II shell evolution \cite{Togashi}. As shown in Fig. \ref{occp} the calculations, clearly indicate that all the relevant states above $0^+_1$ involve on the average three proton excitations from the $1p0f$ orbitals to the $0g_{9/2}$ orbital and five neutron excitations from the $1d_{5/2}2s_{1/2}$ orbitals to the $0g_{7/2}$, $1d_{3/2}$ and $0h_{11/2}$ orbitals. The structure of the $0^+_2$ and $2^+_1$ states are similar to each other. The $0^+_3$ state and all the other states built on top of it also have similar configurations. The $B(E2)$ values within the band predicted from MCSM calculations are very well in agreement with the experiment for the structure based on $0^+_3$.
However, it is to be noted that the MCSM calculations
do not reproduce the strong mixing of the two $2^+$ states (see Fig. \ref{be2_ls}). The observed discrepancies suggest a need for further refinement of the shell model Hamiltonian used in the MCSM calculations.

The five-dimensional beyond mean-field calculations based on the Gogny force (5DCH/D1S) overestimate in general the deformation in $^{98}$Zr.
These calculations suggest a $K=0$ yrast band structure based on $0^+_1$ state, which gets stretched with increasing angular momentum,
i.e. with mean shape coordinates $\beta$ and $\gamma$ changing from 0.23 to 0.37 and 27$^\circ$ to 16$^\circ$, respectively.
The calculations also predict a second well-deformed prolate band including the levels $0^+_2$, $2^+_3$, $4^+_3$, $6^+_3$, $8^+_3$.
This band is predicted to be more rigid against triaxiality ($\gamma\sim18^\circ$); furthermore, the mean axial deformation $\beta$ is changing from 0.34 to 0.42.
The transition strengths within the excited prolate band are at low spin rather close to the experimental values of the band built on the 0$_3^+$ state,
but increase more strongly within the band.
The deformed structures are therefore reasonably well reproduced by this model, however it fails to reproduce the spherical ground state of $^{98}$Zr,
a feature which was also observed in other mass regions and also explains the smooth onset of collectivity at $N=60$ observed for all mean field calculations (see Fig.~\ref{fig4}).

To summarize, using the \RDDS~method~on isotopically identified fission fragments, we have measured the
lifetimes in the neutron-rich isotope $^{98}$Zr, which is located just below the predicted quantum phase transition at $N=60$.
The result shows a very small B($E2;2_1^+\rightarrow0_1^+$) value in $^{98}$Zr confirming the very sudden onset of collectivity at $N=60$.
This effect is well described by the Monte-Carlo Shell-Model and interpreted as Quantum Phase Transition.
Beyond-Mean-Field calculations are in general not able to reproduce this effect, which could be related to the treatment of the particle number projection. The results from the present measurement confirms coexistence of three distinct shapes in this nucleus for the first time. A comparison of the measured $B(E2)$ values with the state-of-art MCSM calculations indicate spherical-prolate-triaxial shape coexistence in $^{98}$Zr.

The authors thank the technical teams at GANIL for their support during the experiment and the authors of Ref. \cite{Navin} for preparing and optimising the VAMOS++ and EXOGAM spectrometers. We also thank A. Navin for helpful discussions. This work has  received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk\l{}odowska-Curie grant agreement No.~702590. It was also supported by the European Community FP7 Integrated Infrastructure Initiative, contract ENSAR No.~262010, by the Research Council of Norway under project grant 213442  and 263030, and by Deutsche Forschungsge-meinschaft (DFG) project DE 1516/3-1. The MCSM calculations were performed on K computer at RIKEN AICS (hp150224, hp160211).
This work was supported in part by the HPCI Strategic Program (The origin of matter and the universe)
and "Priority Issue on post-K computer" (Elucidation of the Fundamental Laws and Evolution of the Universe) from MEXT and JICFuS.