\section{Introduction}

Bayesian inference, which derives the posterior probability using a prior probability and a likelihood function for the observed data, has been an important tool in  statistics.
This approach has been used extensively  for image denoising from early ages to the modern era of deep learning.
For example, in the recent unsupervised deep learning approach for image denoising using Stein's risk
estimate (SURE) \cite{soltanayev2018training}, the unknown Bayesian risk is replaced by the
SURE that can be calculated from the noisy measurement so that deep neural network training is performed by minimizing it.
Unfortunately, this method is sensitive to hyper-parameters, and the neural network must be retrained
if the underlying noise model varies \cite{stein1981estimation}.
On the other hand, there has been increased research interest in  image denoisers that can be trained by minimizing variants of empirical risks that are not associated with clean data.
Noise2Noise ~\cite{lehtinen2018noise2noise} was the first representative approach that does not require clean data.
Unfortunately, multiple noisy versions of the same images are necessary for training.
To address this,  self-supervised learning approaches such as
Noise2Void~\cite{krull2019noise2void}, Noise2Self~\cite{xie2020noise2same}, etc.  have been developed in order to
use only a single  noisy image.
This class of approaches, which  we will call Noise2X throughout the paper, are especially important for practical applications, where noiseless clean images or multiple noisy realization of the same image are difficult or impossible to collect.

Though SURE and Noise2X hold promise for practical applications,
 one of the fundamental questions is how these seemingly different approaches are related and whether there is a coherent mathematical theory that can be leveraged to further improve the performance.
Although we are aware of a preliminary prior work \cite{zhussip2019extending} attempting to unite a subset of them, there is largely a lack of a principled method of designing image denoiser without clean data.

In this article, we try to approach this open problem  from Bayesian statistics in a completely different route.
Instead of minimizing  different forms of empirical risks or SURE,
one of the most important contributions of this paper is the discovery of the importance of classical result from Bayesian statistics  - the Tweedie's formula \cite{efron2011tweedie} which provide an explicit way of computing the posterior mean of canonical parameters from the noisy measurements corrupted with exponential family noises.
Specifically, we show that the Tweedie's formula provides a unified approach for image denoising from any exponential family noises  through
a score function (i.e. the gradient of the loglikelihood).
Therefore, the self-supervised image denoising problem without clean image
can be reduced to the problem of estimating the score function.

\begin{figure}[!t]
	\centering
	\includegraphics[width=0.9\linewidth]{concept.png}

	\caption{ Overall reconstruction flow of (a) supervised learning,  SURE and Noise2X, where  the target $z$ and the distance measure $d(\cdot,\cdot)$
	are uniquely determined by each algorithm, and (b) Noise2Score, where the first step is the estimation of the score function $\hat{l}'$ using neural network training, which is followed by
	Tweedie's formula to obtain the final denoising result.}
	\label{fig:concept}
\end{figure}

This change of perspective has many important theoretical implications as well as the flexibility with regard to the implementation of the algorithm.
 In fact, the score function estimation problem
has been an important research topic in Bayesian statistics and machine learning \cite{hyvarinen2005estimation,vincent2011connection,alain2014regularized}.
In particular,  Alain and Bengio~\cite{alain2014regularized} showed that the minimization of the denoising autoencoder (DAE) objective function
provides  an explicit way of approximating  the score function. This result was further
 extended using the amortized residual denoising autoencoder  (AR-DAE) for numerical stability and accuracy \cite{lim2020ar}.
 Interestingly, the training method for DAE or its variants is strikingly similar to that of Noise2X. Therefore,  by combining this with the Tweedie's formula, we can obtain a novel unified
 framework of Bayesian approach for self-supervised image denoising.
The conceptual difference and similarity of our method to the existing approaches are illustrated in Fig.~\ref{fig:concept}(a)(b).
As the first step of our algorithm   is similar to Noise2X, we call our method Noise2Score.

Noise2Score is so powerful that it can be used to deal with any exponential family noises.
Moreover,  in contrast to the SURE approach, in which the network must be retrained using a different loss function depending on the noise model and parameters \cite{soltanayev2018training,kim2020unsupervised}, Noise2Score is universal in the sense that the identical neural network training step is used regardless of noise models.
 This property can be exploited to extend Noise2Score to blind setups in which noise parameters are unknown and should be estimated with minimal complexity.
In addition to the novel theoretical findings mentioned above, our empirical results using additive Gaussian, Poisson and Gamma noise models have shown that Noise2Score significantly outperforms SURE and Noise2X under similar experimental conditions.

\section{Related Works}

Here we give a brief overview of the existing deep learning approaches to denoise images.
As shown below, their main differences arise from the different choices of the distance metric $d(\cdot, \cdot)$
in Fig.~\ref{fig:concept}.

\subsection{Image denoising using supervised learning}

Let $\mathcal{X}$ and $\mathcal{Y}$ denotes the spaces of clean and noisy images with the probability
measure $P_X$ and $P_Y$, respectively.
When the paired ground-truth data $x$ is available for a noisy measurement $y$, i.e. $x:=x(y)$, then
 the supervised training can be carried out by minimizing the following loss function:
 \begin{eqnarray}
\ell_{sup}(\Theta) = {{\mathbb E}}_{y\sim P_Y}d(x, F_\Theta(y)),&\mbox{where} & d(x,F_\Theta(y)):=\|x - F_\Theta(y)\|^2
\end{eqnarray}
where $F_\Theta(y)$ is a deep neural network parameterized by $\Theta$ with input $y$ and
here we use $l_2$ distance metric to simplify the explanation.
Although this approach has significantly improved performance compared to the classical denoising algorithms \cite{dabov2006image},
the supervised learning approach has a drawback as the large number of paired clean images are necessary.
In real world applications, the acquisition for paired clean images are difficult or even impossible.

\subsection{Noise2X for image denoising without clean images}

To address this issue, various forms of the denoisers that can be trained without clean images have been proposed, which includes
Noise2Noise~\cite{lehtinen2018noise2noise}, Noise2Void~\cite{krull2019noise2void}, Noise2Self~\cite{batson2019noise2self}, Noise2Same ~\cite{xie2020noise2same}, etc.
These approaches are usually based on the variants of  loss functions that do not depend on clean images:
\begin{eqnarray}
\ell_{Noise2X}(\Theta) = {{\mathbb E}}_{y\sim P_Y}d(y', F_\Theta(y)),&\mbox{where} & d(y',F_\Theta(y)):=\|y' - F_\Theta(y)\|^2
\end{eqnarray}
where the target image $y'$ is related to $y$ in unique ways depending on algorithms.
For example,  in Noise2Noise~\cite{lehtinen2018noise2noise}, $y'$ is
 another noisy realization of the same underlying clean images. However,  Noise2Noise is not practical
 since multiple noisy realizations of the same image are not usually available in practice.
 Noise2Void~\cite{batson2019noise2self} generates $y'$ by adding blind spots on $y$.
 In Noise2Self ~\cite{batson2019noise2self}, so-called J-invariant mask $J$  is added on the input $y$,
and $y'$ is defined as the image on the complementary mask.
All of these changes of the target images are intended to prevent the network from converging to a trivial identity mapping
\subsection{Image denoising using Stein's Unbiased Risk Estimate (SURE)}

Rather than using the variants of a  loss function as in Noise2X,
Soltanayev et al~\cite{soltanayev2018training} proposed a training scheme using the loss function from
Stein's unbiased risk estimate (SURE) \cite{stein1981estimation}
which has additional regularization term.
Specifically, the loss function for the SURE denoiser is defined by
\begin{eqnarray}\label{eq:SURE}
\ell_{SURE}(\Theta) =  {{\mathbb E}}_{y\sim P_Y} d(y, F_\Theta(y)),&~ d(y,F_\Theta(y)):= \|y- F_\Theta(y)\|^2+ 2\sigma^2 \mathrm{div}_y F_\Theta(y)
\end{eqnarray}
Here,  the second term in $d(y,F_\Theta(y))$ of \eqref{eq:SURE} is the regularization term to prevent the network from converging to a trivial identity mapping,
where $\mathrm{div}_y$ denotes the divergence with respect to $y$.
Due to the difficulty of calculating the divergence term, the authors in \cite{soltanayev2018training} relied on MonteCarlo SURE \citep{ramani2008monte} which calculates the divergence term using MonteCarlo simulation.
This introduces additional hyperparameters, on which the final results critically depend.
Yet another limitation of SURE image denoiser is that the loss function is determined by the noise model.
For example,
in \cite{kim2020unsupervised}, the authors derived Poisson Unbiased risk estimator (PURE) to train deep CNN for dealing with Poisson noise,
whose loss function is completely different from \eqref{eq:SURE}. Even under the same Gaussian noise models, if the noise variance $\sigma^2$ varies,
the neural network $F_\Theta$ must be trained again, which causes additional computational burden.
This differs from Noise2X, which uses the same loss function regardless of the noise models.

\section{Main Contribution: Noise2Score}

In contrast to the existing approaches shown in Fig.~\ref{fig:concept}(a),  our denoising approach  is composed of two steps  as illustrated in Fig.~\ref{fig:concept}(b),  which is inspired by the results from Bayesian statistics.
In the following, we explain the details of each step.

\subsection{Tweedie's formula of the posterior mean for image denoising}

Suppose that the noisy measurement $y$ is given by
\begin{eqnarray}\label{eq:gaussian}
y=x+w,&\mbox{where}& w \sim \mathcal{N}(0,\sigma^2I)
\end{eqnarray}
where $x$ is the clean image, and $\mathcal{N}(0,\sigma^2I)$ denotes
the independent and identically distributed  (i.i.d) Gaussian distribution with zero mean and variance $\sigma^2$.
Then, Tweedie’s formula calculates the posterior expectation of $x$ given $y$ as \cite{robbins1956empirical}
\begin{eqnarray}\label{eq:post}
{{\mathbb E}}[x|y]= y+\sigma^2 l'(y),&\mbox{where}& l'(y)=\nabla_y \log p(y)
\end{eqnarray}
where $p(y)$ denotes the marginal distribution of $y$.
For the image corrupted by   Gaussian noises,
the posterior mean is the  minimum mean square error (MMSE) estimator, so that
 the denoised image can
be obtained  using \eqref{eq:post}, as long as
we know $l'(y)$ - the {\em score function} \cite{serfling2009approximation}.

Tweedie's formula was extended for general class of exponential family distribution \cite{efron2011tweedie}.
Specifically,
exponential family of probability distributions are defined as those
distributions whose density  have the following general form:
\begin{eqnarray}
p(y|\eta)=p_0(y)\exp\left(\eta^\top T(y) -\varphi(\eta) \right)
\end{eqnarray}
where the superscript $^\top$ denotes the transpose operation.
Here $\eta$ is a canonical (vector) parameter  of the family, $T(y)$ is a (vector) function of $y$,
 $\varphi(\eta)$ the cumulant generating
function which makes $p(y|\eta)$
 integrate to 1, and $p_0(y)$ the density up to a scale factor when $\eta = 0$.
Exponential family distribution covers a large class of important distributions such as
the Gaussian, binomial, multinomial,
Poisson, gamma, and beta distributions, as well as many others.
Table~\ref{tbl:exp} summarizes the corresponding parameters for each exponential family distribution.

\begin{table*}[!b]
	\centering
	\caption{Tweedie's formula of exponential family distribution for image denoising. }
	\label{tbl:exp}
	\resizebox{\textwidth}{!}{
		\begin{tabular}{c|ccccc|c}
			\hline
			{Distribution}					& $p(y|x)$			&  $\eta$ & $T(y)$	 & $p_0(y)$   & $l_0'(y)$ & $\hat x$	\\ \hline\hline
			Gaussian  &  		$\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(y-x)^2}{2\sigma^2}}$				& ${x}/{\sigma^2}$	& $y$	 &$\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{y^2}{2\sigma^2}}$	 	 & $-\frac{y}{\sigma^2}$  & $y+\sigma^2 l'(y)$ \\
			Gaussian  &  		$\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(y-x)^2}{2\sigma^2}}$				& $\left[\frac{x}{\sigma^2}, -\frac{1}{2\sigma^2} \right]^\top$	& $\left[y,y^2\right]^\top$	 &$\frac{1}{\sqrt{2\pi}} $	 	 & $0$  & $y+\sigma^2 l'(y)$ \\
			Poisson	&	$\frac{x^ye^{-x}}{y!}$	& $\log(x)$	& $y$ 	&	  $\frac{1}{y!}$  & $\simeq -\log\left(y+\frac{1}{2}\right)$ &$\left(y+\frac{1}{2}\right)\exp(l'(y))$ \\
			Gamma($\alpha,\beta$) & $\frac{\beta^\alpha}{\Gamma(\alpha)}\left(\frac{y}{x}\right)^{\alpha-1}e^{-\beta\frac{y}{x}}$  & $\left[\alpha-1,-\frac{\beta}{x}\right]^\top$ & $\left[\log y,-y\right]^\top$ &  1  &0
			& $\frac{\beta y}{(\alpha-1)-y l'(y)}$ \\
			Bernoulli & $x^y(1-x)^{(1-y)}$ & $\log\left(\frac{x}{1-x}\right)$ &  $y$ & $1$ & 0 & $\frac{e^{l'(y)}}{1+e^{l'(y)}}$ \\
			Exponential & $xe^{-yx},~ y\geq 0$ &  $-x$ & $y$    & 1  & 0  & $-l'(y)$\\

			\hline
		\end{tabular}
	}
\end{table*}

Then, using the Bayes' rule, the posterior density of $\eta$ given $y$ is given by \cite{efron2011tweedie}:
\begin{eqnarray}
p(\eta|y) = \exp\left(\eta^\top T(y) -\lambda(y) \right)\left[p(\eta)e^{-\varphi(\eta)}\right],&\mbox{where}& \lambda(y)=\log\left(\frac{p(y)}{p_0(y)}\right)
\end{eqnarray}
where $p(y)$ and $p(\eta)$ denote the marginal distribution of $y$ and $\eta$, respectively.
This implies that the posterior density is again an exponential family distribution,
and the mode of the  posterior distribution  can be obtained by finding the maximum of $p(\eta|y)$.
Specifically, by computing the gradient of $\log p(\eta|y)$ with respect to $y$ and setting it to zero,
the posterior estimate of the canonical parameter $\hat \eta$ should satisfy the following equality:
\begin{eqnarray}\label{eq:fix}
\hat\eta^\top  T'(y) &=&\lambda'(y)=-\nabla_y \log p_0(y)+\nabla_y \log p(y)
=  -l_0'(y)+l'(y)
\end{eqnarray}
where $l'(y):=\nabla_y \log p(y)$ and $l_0'(y):=\nabla_y \log p_0(y)$ are score functions, and $T'(y)=\nabla_y T(y)$.
In particular, if $T(y)=y$, we can obtain the following closed form solution for the posterior mean:
\begin{eqnarray}
\hat\eta:={{\mathbb E}}[\eta|y]= \lambda'(y
& =&  -l_0'(y)+l'(y)
\end{eqnarray}
In the following, we derive denoising algorithms using Tweedie's formulation  for several representative
distributions.

\paragraph{Additive Gaussian noise}

For  Gaussian distribution, in the first row of Table~\ref{tbl:exp}, we have
\begin{eqnarray}
p_0(y)=\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{y^2}{2\sigma^2}}, && \eta = \frac{x}{\sigma^2}
\end{eqnarray}
Accordingly,
\begin{eqnarray}
{{\mathbb E}}[\eta|y]= \frac{\hat x}{\sigma^2} =  \frac{y}{\sigma^2} + l'(y), &\mbox{where} & \hat x:= {{\mathbb E}}[x|y]
\end{eqnarray}
since $ l_0'(y)={y}/\sigma^2$. Therefore, the posterior mean  is given by
\begin{eqnarray}
\hat x = y + \sigma^2 l'(y)
\end{eqnarray}
which is equal to \eqref{eq:post}. Even with the different parameterization of Gaussian as shown
in the second row of Table~\ref{tbl:exp}, we can arrive at the same result using \eqref{eq:fix}.

\paragraph{Poisson noise}

For the case of Poisson noises, we have
\begin{eqnarray}
p_0(y)=\frac{1}{y!}&, &l_0'(y) =  -\frac{\nabla_y \Gamma(y+1)}{\Gamma(y+1)} \simeq -\log\left(y+\frac{1}{2}\right)
\end{eqnarray}
where $\Gamma(y)$ denotes the gamma function and the last approximation comes from \cite{efron2011tweedie}.
Accordingly, we have
\begin{eqnarray}
\hat \eta = \log (\hat x)=  \log\left(y+\frac{1}{2}\right)  + l'(y)
\end{eqnarray}
This leads to  the following posterior estimate of the image
\begin{eqnarray}
\hat x = \left(y+\frac{1}{2}\right)\exp(l'(y))
\end{eqnarray}
This corresponds to Tweedie's formula for Poisson case~\cite{robbins1956empirical}.

This result can be generalized to the low photon count sensing scenario, where the sensor
measurement can be described by \cite{le2014unbiased}:
\begin{equation}
    y = \zeta z, \; z \sim \text{Poisson}(x/\zeta)
\label{poisson}
\end{equation}
where $z$ is a random variable that follows the Poisson distribution in Table~\ref{tbl:exp}
 and $\zeta >0$ is the gain of the acquisition process that is related to the noise level.
 Then, we have
 \begin{eqnarray}
p_0(z)=\frac{1}{z!}, && \eta = \log\left(\frac{x}{\zeta}\right)
\end{eqnarray}
Therefore, we have
\begin{eqnarray}
\frac{\hat x}{\zeta} =\left(z+\frac{1}{2}\right)\exp(l'(z)) & \longrightarrow &  {\hat x} =\left(y+\frac{\zeta}{2}\right)e^{ l'\left(\frac{y}{\zeta}\right)}
\end{eqnarray}

\paragraph{Gamma noise}

Gamma noise distribution can be used to model the speckle noises in various imaging application.
Specifically, the image corrupted with the speckle noise is represented by
\begin{equation}
y = x n, \quad  n\sim p(n;\alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}n^{\alpha-1}\exp(-\beta n)
\end{equation}
where $p(n;\alpha,\beta)$ denotes the Gamma distribution with $(\alpha,\beta)$ parameters.
Here,  $(\alpha,\beta)$ are the parameters which determine the noise level of Gamma noise distribution.
For example, for the $k$-look measurement case, $\alpha=\beta=k$ \cite{bioucas2010multiplicative}.
Then, the probability $p(y|x)$ can be obtained by replacing $n$ in $p(n;\alpha,\beta)$ with
$n=y/x$, and the resulting density
function can be bound in Table~\ref{tbl:exp}.
Then, we can easily show that
\begin{equation}
  p_o(y) = 1, \quad
 \eta =
\begin{bmatrix}
\alpha-1 \\
-\beta/x
\end{bmatrix}, \quad
T(y) =
\begin{bmatrix}
\log y\\
y
\end{bmatrix}
\label{gamma}
\end{equation}
Using \eqref{eq:fix}, we have
\begin{eqnarray}\label{eq:gamma}
\frac{\alpha-1}{y} -\frac{\beta}{x} = l'(y)& \Longrightarrow &
\hat x = \frac{\beta y}{(\alpha-1)-yl'(y)}
\end{eqnarray}

Similar derivations for other exponential family distribution
can be found in Table~\ref{tbl:exp}.

\subsection{Score function estimation}

So far, our derivation assumes the prior knowledge of the score function.
In  practice, this should be estimated, so here we describe how the score function can be estimated using a neural network.

Historically, Hyv\"{a}rinen et al
\cite{hyvarinen2005estimation}  was the first to derive a remarkable implicit
score matching objective  that no longer requires having an explicit score target  but is nevertheless equivalent to the original
problem that minimizes the expected quadratic distance between the model function and the score function of data. With the advance of the Denoising Auto Encoder (DAE) \cite{vincent2010stacked}, the author of \cite{vincent2011connection} observed that the minimization of the DAE objective function with the residual form is related to the score matching between the model and perturbed data.
This observation was rigorously analyzed by
 Alain and Bengio~\cite{alain2014regularized}, who showed that
as the perturbed noise becomes sufficiently small, DAE leads to the score function

More specifically, in the denoising autoencoder (DAE), the following loss is minimized:
\begin{equation}\label{eq:DAEcost}
    \ell_{DAE}(\Theta) =
    \underset{\sigma_a \sim \mathcal{N}(0,\delta^2)}{\underset{u \sim \mathcal{N}(0,I)}{\underset{y \sim P_Y}{{{\mathbb E}}}}} \|y -F_\Theta(y+\sigma_a u)\|^2
\end{equation}
According to  \cite{alain2014regularized}, the optimal DAE $F_{\Theta^*}(x)$
can be represented by
\begin{equation}
F_{\Theta^*}(y) = y + \sigma_a^2 l'(y) + o(\sigma_a^2),
\label{bengio}
\end{equation}
where $l'(y)$ is the score function defined in \eqref{eq:post}
and $o(\cdot)$ denotes the small ``o'' notation.
In other words, for sufficiently small $\sigma_a$, we can approximate the score function as:
\begin{eqnarray}\label{eq:DAEs}
\hat l'(y)= \frac{F_{\Theta^*}(y)-y}{\sigma_a^2}
\end{eqnarray}

One of the downsides of the score function estimation using DAE is that
the score function estimate in \eqref{eq:DAEs} is numerically unstable as $\sigma_a \rightarrow 0$.
To address the numerical instability and reduce the approximation error,
Lim et al \cite{lim2020ar} recently proposed so-called the amortized residual DAE (AR-DAE).
Specifically, AR-DAE is trained  by minimizing the following the objective function:
\begin{equation}
    \ell_{AR-DAE}(\Theta) = \underset{\sigma_a \sim \mathcal{N}(0,\delta^2)}{\underset{u \sim \mathcal{N}(0,I)}{\underset{y \sim P_Y}{{{\mathbb E}}}}}\|u + \sigma_a R_\Theta(y + \sigma_a u)\|^2
\label{eq:ardae}
\end{equation}
where
$R_\Theta$ is  from the residual form of the DAE:
\begin{eqnarray}\label{eq:res}
F_\Theta(y)=\sigma_a^2 R_\Theta(y)+y
\end{eqnarray}
By plugging  \eqref{eq:res} in \eqref{eq:ardae}, we can obtain the original DAE cost in \eqref{eq:DAEcost} up to a scale factor.
Furthermore, using \eqref{eq:DAEs}, we have
\begin{eqnarray}\label{eq:gaus}
\hat l'(y)= \frac{F_{\Theta^*}(y)-y}{\sigma_a^2} = R_\Theta(y)
\end{eqnarray}
Therefore,
the neural network trained with the AR-DAE
\eqref{eq:ardae} is a direct and stable way of estimating the score function.
Therefore, we employ the AR-DAE to estimate the score function from the noisy measurement as the first step of our method
(see Fig.~\ref{fig:concept}(b)).

\subsection{Relation to Noise2X and SURE}

Although our Noise2Score was derived from a completely different perspective, it turns out that it has very important connections to Noise2X and SURE.

Specifically, if the noisy image $y$ is corrupted by additive Gaussian,
using Tweedie's fomula for Gaussian noises in Table~\ref{tbl:exp}, the denoised image can be approximated by
\begin{eqnarray}\label{eq:xhat}
\hat x = y + \sigma^2 \hat l'(y) = F_{\Theta^*}(y)
\end{eqnarray}
where we use \eqref{eq:gaus}.
Recall that $F_{\Theta^*}(y)$ is the neural nework trained with DAE loss function in \eqref{eq:DAEcost}.
The DAE training with \eqref{eq:DAEcost} is basically inserting noises into the images and then find the mapping that removes the noises.
\add{In terms of adding extra noise to a noisy image, our method is closely related to Noisier2Noise~\cite{moran2020noisier2noise}. However, in the training phase of Noisier2Noise, the main assumption is that noise from the same noise statistic as the original noisy image should be injected into the noisy images so that the neural networks learn to reduce the noise from the noisy input image. In our Noise2Score training method, however, Gaussian noise is added to noisy images at different noise levels regardless of the noise statistics to estimate the score function rather than noise. In addition, in contrast to Noisier2Noise, our method has a post-processing step with the Tweedie’s formula, which is determined by the noise model such as Poisson, Gamma, etc.}

Yet another important connection is its relation to SURE.
Using the residual representation in \eqref{eq:res}, the
SURE cost function in \eqref{eq:SURE} can be equivalently
represented by
\begin{eqnarray}
\ell_{SURE}(\Theta) &= &  {{\mathbb E}}_{y\sim P_Y} \left\{\|y- F_\Theta(y)\|^2+ 2\sigma^2 \mathrm{div}_y F_\Theta(y)  \right\} \notag \\
&= & {{\mathbb E}}_{y\sim P_Y} \left\{\sigma^4 \|R_\Theta(y)\|^2+ 2\sigma^4 \mathrm{div}_y R_\Theta(y)  \right\} +2\sigma^2 \mathrm{dim}(y)\
\end{eqnarray}
where $\mathrm{dim}(y)$ is the dimension of the vector $y$, which is constant.
Up to the scaling factor, this cost function is identical to the
implicit score matching objective by Hyv\"{a}rinen et al
\cite{hyvarinen2005estimation}, which is given by
\begin{eqnarray}\label{eq:ISM}
\ell_{ISM}(\Theta) =  {{\mathbb E}}_{y\sim P_Y} \left\{ \frac{1}{2} \|\Psi_\Theta(y)\|^2+ \mathrm{div}_y \Psi_\Theta(y) \right\}
\end{eqnarray}
where $\Psi_\Theta(y)$ is the score function estimate parameterized by $\Theta$.
Therefore, using the residual transform in \eqref{eq:res}, the SURE objective is nothing but to find
the score function, which corresponds to the first step of our Noise2Score. Furthermore, due to the use of \eqref{eq:res}, SURE is not optimal for noise models other than Gaussian. \add{Although there are variants of SURE for other noise type~\cite{raphan2011least}, the cost function of this method varies depending on the type of noise. However, as discussed in detail later, our method uses the same cost function for neural network training, but only the post-processing step differs depending on the type of noise.}
Moreover, considering that the modern research trend of the score function estimation has quickly evolved from the implicit score matching
to DAE or AR-DAE due to the numerical instability
and difficulty of computing the divergence term in \eqref{eq:ISM}, we can easily expect that
the SURE denoiser also suffers from similar limitations.

\subsection{Universal neural estimation for blind denosing problems}

Recall that Noise2Score is composed of two steps: the score function estimation by minimizing the AR-DAE loss function,
after which the final denoised image is obtained by applying Tweedie's formula  shown in Table~\ref{tbl:exp}.
This decoupling allows Noise2Score to have unique advantages compared to the conventional  approaches.

Specifically, one of the most important advantages of Noise2Score is that
AR-DAE training is universal in the sense that
the same loss function is used regardless of noise models and parameters.
More specifically,  in  \eqref{eq:ardae},  the model
and noise parameter dependency is only through the sampling $y$ from
 $P_Y$. This sampling  step can be performed by random patch cropping from the input images, which procedure
 is identical regardless of $P_Y$ being Gaussian,  Poisson,  Gamma, etc.

Yet another advantage is that this property can be utilized to estimate the unknown noise parameters  without retraining neural networks.
For example, in the case of Gaussian denoising in Table~\ref{tbl:exp},
we have the  following form of the denoised images:
\begin{eqnarray}
\hat x(\sigma)= y +\sigma^2 \hat l'(y)
\end{eqnarray}
If the noise parameter $\sigma$ is unknown, then we can
find the optimal $\sigma^*$ value by solving the optimization problem
\begin{eqnarray}\label{eq:Q}
\sigma^*=\arg\min\limits_\sigma Q\left(\hat x(\sigma)\right)
\end{eqnarray}
where $Q(\cdot)$ is an image quality penalty such as total variation (TV).
In contrast to the SURE, where the $\sigma$ dependency exists in the loss function in \eqref{eq:SURE} and
the neural network must be trained with a new $\sigma$ value,
the optimization problem in \eqref{eq:Q} is  computationally negligible since the precalculated score function estimate $\hat l'(y)$ is
used regardless of different $\sigma$ values.
A similar on-the-fly parameter estimation can be used for Gamma noise model. As shown in \eqref{eq:gamma}, the parameters
$(\alpha,\beta)$ can be adjusted using the precomputed score function estimate $\hat l'(y)$. Even for the Poisson noise model, the same trained neural network is used to compute the score function for each input $y/\zeta$.
\add{Therefore, our Noise2Score has advantages for the blind application where  existing methods  cannot be used due to the difficulty of  generating samples by changing the noise level parameters without retraining.}

\section{Experimental Results}
\paragraph{Dataset and Implementation detail}
We evaluated the proposed method for gray-scale and color images in the four benchmark datasets: the gray-scale image dataset contains BSD68~\cite{martin2001database} and Set12. The color image dataset contains RGB natural images CBSD68~\cite{martin2001database}, Kodak dataset. We adopt DIV2K~\cite{Timofte_2018_CVPR_Workshops}, and (C)BSD400 dataset as train data set. We generated the synthetic noise images for each noise distribution. In order to evaluate the proposed method fairly with comparison methods,  the same modified U-Net generator \cite{lehtinen2018noise2noise} is used for all methods. The total epoch was set to 100 and the Adam optimizer~\cite{kingma2014adam} was used to train the network. The learning rate was initialized to 2$\times10^{-4}$ for first 50 epoch, and after 50 epoch the learning rate was decayed to $2\times10^{-5}$. The proposed method was implemented in PyTorch~\cite{paszke2017automatic} with NVidia GeForce GTX 1080-Ti. The network training  took about 10 hours. \add{The more detail of implementation are described in Supplementary Material}. {To deal with blind noise case,  as for
the quality metric $Q(\cdot)$ we use the TV norm for the Gaussian noises and its combination with  data fidelity term for the Poisson and Gamma noise. The
details of the quality metric and the discussion on the accuracy of the blind parameter estimation is provided in  Supplementary Material.}

\begin{table}[t]
	\begin{small}
\begin{center}
			\caption{Quantitative comparison for various noise model using various methods in terms of PSNR(dB) when the noise parameters are known or unknown
			(N2V: Noise2Void, N2S: Noise2Self, N2N: Noise2Noise, SL: supervised learning, Anscombe: Anscombe method for BM3D).}
				\label{tbl-results}
\resizebox{\textwidth}{!}
{%
\begin{tabular}{clcccccccccc}\toprule
	\multicolumn{1}{c}{Noise type}                    & & \multicolumn{7}{c}{Known parameters}                                     & \multicolumn{3}{c}{Unknown  parameters} \\ \cmidrule(r){3-9} \cmidrule(r){10-12}
	Gaussian                    & Dataset & BM3D     & N2V   & N2S   & SURE  & Ours           & N2N   & SL   & Ours        & N2N & SL   \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$\sigma=$ 25}  & BSD68   & 28.59    & 26.77 & 28.28 & 28.78 & \textbf{29.12} & 29.18 & 29.20 & 28.95       & 28.98  & 28.98      \\
	& Set12   & 29.96    & 27.56 & 29.16 & 29.13 & \textbf{30.13} & 30.33 & 30.36 & 30.08       & 30.08   & 30.08    \\
	& CBSD68  & 30.56    & 29.22 & 30.05 & 30.23 & \textbf{30.85} & 31.10 & 31.10 & 30.78       & 30.91   & 30.91    \\
	& Kodak   & 31.68    & 30.02 & 30.53 & 30.75 & \textbf{31.89} & 32.20 & 32.20 & 31.78       & 31.96   & 31.96     \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$\sigma$ = 50}  & BSD68   & 25.62    & 24.34 & 25.61 & 25.80 & \textbf{26.21} & 26.27 & 26.30 & 25.81       & 25.86  & 25.98     \\
	& Set12   & 26.33    & 24.68 & 26.19 & 26.23 & \textbf{27.16} & 27.20 & 27.20 & 26.65       & 26.65  & 26.72    \\
	& CBSD68  & 27.38    & 25.13 & 27.05 & 26.24 & \textbf{27.75} & 27.94 & 27.95 & 27.32       & 27.66  & 27.68     \\
	& Kodak   & 27.02    & 25.75 & 28.01 & 26.93 & \textbf{28.83} & 29.07 & 29.10 & 28.13       & 28.70  & 28.71     \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	Poisson                    & Dataset & Anscombe & N2V   & N2S   & PURE  & Ours           & N2N   & SL   & Ours  & N2N   & SL         \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$\zeta$ = 0.01} & BSD68   & 30.51    & 28.73 & 29.76 & 30.16 & \textbf{30.81} & 30.89 & 30.98 & 30.63       & 30.91  & 30.93     \\
	& Set12   & 31.20    & 30.06 & 30.47 & 30.69 & \textbf{31.58} & 31.70 & 31.79 & 31.42       & 31.67  & 31.70     \\
	& CBSD68  & 32.40    & 31.85 & 31.04 & 32.30 & \textbf{32.61} & 33.01 & 33.01 & 32.23       & 32.94  & 32.95     \\
	& Kodak   & 33.13    & 32.98 & 32.24 & 33.01 & \textbf{33.41} & 33.91 & 33.91 & 32.96       & 33.85  & 33.87     \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$\zeta$ = 0.05} & BSD68   & 26.77    & 26.12 & 26.54 & 24.76 & \textbf{27.12} & 27.23 & 27.25 & 26.77       & 27.19  & 27.20      \\
	& Set12   & 27.54    & 27.21 & 27.53 & 25.07 & \textbf{27.86} & 28.02 & 28.04 & 27.76       & 28.01  & 28.03     \\
	& CBSD68  & 28.33    & 28.37 & 28.32 & 26.68 & \textbf{28.68} & 29.23 & 29.27 & 28.23       & 29.07  & 29.07     \\
	& Kodak   & 29.31    & 29.56 & 29.52 & 26.38 & \textbf{29.71} & 30.33 & 30.40 & 28.98       & 30.13  & 30.25     \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	Gamma                   & Dataset &       & N2V   & N2S   &     & Ours           & N2N   & SL   & Ours  & N2N  & SL         \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$k$ = 100}    & BSD68   & -        & 29.32 & 30.49 & -     & \textbf{32.67} & 32.87 & 32.93 & 32.54       & 32.83 & 32.88      \\
	& Set12   & -        & 30.54 & 30.71 & -     & \textbf{33.01} & 33.21 & 31.72 & 32.89       & 33.19 & 33.21       \\
	& CBSD68  & -        & 31.11 & 30.54 & -     & \textbf{33.82} & 35.45 & 35.53 & 33.50       & 35.33  & 35.33     \\
	& Kodak   & -        & 31.96 & 31.60 & -     & \textbf{34.22} & 36.26 & 36.41 & 33.82       & 36.16   & 36.16    \\ \cmidrule(r){1-2} \cmidrule(r){3-9} \cmidrule(r){10-12}
	\multirow{4}{*}{$k$ = 50}     & BSD68   & -        & 26.98 & 29.25 & -     & \textbf{30.53} & 31.10 & 31.16 & 30.42       & 30.64 & 30.64      \\
	& Set12   & -        & 27.36 & 29.67 & -     & \textbf{30.87} & 31.68 & 31.72 & 30.83       & 31.58  & 31.67     \\
	& CBSD68  & -        & 30.51 & 30.19 & -     & \textbf{31.05} & 33.52 & 33.53 & 30.93       & 33.40  & 33.41     \\
	& Kodak   & -        & 31.38 & 31.03 & -     & \textbf{31.34} & 34.49 & 34.57 & 31.32       & 34.39  & 34.40   \\ \bottomrule
\end{tabular}}
	\end{center}
	\end{small}
\end{table}

\begin{figure}
    \centering
	\includegraphics[width=0.9\linewidth]{results.jpg}
    \caption{Qualitative comparison using CBSD68 dataset. Top : Gaussian noise  with $\sigma$ = 50. Middle: Poisson noise with $\zeta$ = 0.05. Bottom: Gamma noise with $k$ = 100. White numbers at the lower left part of the images indicate the PSNR values in dB.  Noisy: noisy input,  GT:  ground-truth image.}
        \label{fig:result}
\end{figure}

\paragraph{Gaussian noise}
Using additive Gaussian noise, our method was compared to (C)BM3D \cite{dabov2006image}, Noise2Void, Noise2Self, SURE, Noise2Noise and supervised  learning  as shown in Table 2.
As expected, the supervised learning and  Noise2Noise using multiple noisy realization were the best, but these methods are not practical as explained before.
Among the other self-supervised learning approaches, our Noise2Score provide the best performance for all dataset.
With the comparison with the supervised learning and Noise2Noise, the performance difference of our method is only margin.
Even for blind-noise cases where the noise levels are randomly sampled from Gaussian distribution with $\sigma \in$ [5, 55],
Noise2Score with the on-the-fly parameter estimation
performed comparable to the supervised learning approach. These results indicate the superior performance of our method.
The qualitative comparison in Fig.~\ref{fig:result} show that Noise2Score retains the better image details and provides
much visually pleasing results compared to other top two self-supervised learning approaches.

\paragraph{Poisson noise}
For the Poisson noise case, BM3D and SURE were replaced with Anscombe (BM3D+VST)~\cite{makitalo2010optimal} and  PURE \cite{kim2020unsupervised}, respectively,  as indicated in Table 2, as  they are optimized for Possion noise model.
Again, our Noise2Score results in significant performance gain compared to other approaches,
and the results are comparable to supervised learning approach and Noise2Noise.
In addition, to perform experiments when the noise parameters are unknown, we performed experiments with  the noise levels randomly sampled between $\zeta \in$ [0.001,0.1]. We found that our proposed method produce even higher PSNR than the results of  other existing methods even with known parameters,
and also provide comparable results to the supervised learning approaches.
The qualitative comparison in Fig.~\ref{fig:result} with the other two best approach showed that
our proposed method provides much cleaner images.

\paragraph{Gamma noise}
As the extension of BM3D and SURE are not available for Gamma noises,
only four comparison methods was adopted for comparison as indicated in Table~\ref{tbl-results}.  Here, we set $\alpha=\beta=k$.
Again,
 Noise2Score produced best results among self-supervised learning approaches, and provided  comparable results to  Noise2Noise and supervised learning approaches.
Even for the blind noise case, where the noise parameter $k$ randomly selected in $k \in [40,120]$,  Noise2Score provided the better PSNR results
compared to the other self-supervised learning method even with known parameters, and produces comparable results to supervised learning approaches.
 The qualitative comparison in Fig.~\ref{fig:result} confirm that our Noise2Score provides best reconstruction results.

\section{Conclusion}

In this work, we  provided a novel Bayesian framework for self-supervised im\label{key}age denoising without clean data, which surpasses SURE, PURE, Noise2X, etc.  Our novel innovation came from  the Tweedie's formula,
which provides explicit representation of denoise images through the score function. By combining with the score-function estimation using AR-DAE, our Noise2Score
can be applied to image denoising problem from any exponential family noises. Furthermore,  an identical neural network training can be
universally used regardless of the noise models, which leads to  the noise parameter estimation with minimal complexity.
The links to SURE and existing Noise2X were also explained, which clearly showed why our method is a better generalization.

\section*{Limitation and negative societal impacts}
As a negative societal impact, the failure of image denoising methods could produce side effects. For example, removing both the noise and the texture of the medical images could lead to misdiagnoisis. While the proposed Noise2Score has merits, there are also limitations. In real environments, the prior knowledge of noise distribution may not be available, and the noise model could not be modeled by exponential family noises. Therefore, future extension to such scenario would be beneficial.

\section*{Acknowledgment}
This research was funded by the National Research Foundation (NRF) of Korea grant NRF-2020R1A2B5B03001980,  ETRI (Electronics and Telecommunications Research Institute)’s internal funds [21YR2500, Development of Digital Biopsy Core Technology for high-precision Diagnosis and Therapy of Senile Disease], and KAIST Key Research Institute (Interdisciplinary Research Group) Project.
This work was also supported by Institute of Information $\&$ Communications Technology Planning $\&$ Evaluation (IITP) grant funded by the Korea government (MSIT) (No.2019-0-00075, Artificial Intelligence Graduate School Program (KAIST)).

{\small
\bibliographystyle{unsrt}


\section*{Introduction}
\label{Introduction}

Machine learning techniques have demonstrated remarkable success
in numerous topics in science and engineering, including artificial intelligence \cite{mitchell1990machine,duda2012pattern},  molecular dynamics \cite{botu2015adaptive}, light harvesting systems \cite{hase2017machine}, molecular electronic properties \cite{montavon2013machine}, surface reaction network \cite{ulissi2017address}, density functional models \cite{brockherde2017bypassing}, phase classification, and quantum simulations \cite{wang2016discovering,carrasquilla2017machine,broecker2017machine,ch2017machine,van2017learning,arsenault2014machine,kusne2014fly}. In addition, modern machine learning techniques have also been applied to the state space of complex condensed-matter systems for their abilities to analyze exponentially large data sets \cite{carrasquilla2017machine},  speed-up searches for novel energy generation/storage materials \cite{de2017use,wei2016neural} and classification of  entanglement\cite{gao2018experimental}. With the rapid development of quantum computers\cite{leibfried2003experimental,debnath2016demonstration, karra2016prospects, arute2019quantum, zhong2020quantum}, it has become a new frontier to recognize patterns using quantum computers. Considering recent advancements in both quantum computing and machine learning, the combination of the two techniques – quantum machine learning – is expected to be a promising application of quantum computer in the near future. Many quantum machine learning algorithms were proposed in the past few years\cite{rebentrost2018quantum,cao2016combined,biamonte2017quantum,roy2021enhancement,dixit2021training}.
Moreover, researchers have succeeded to apply quantum machine learning  algorithms to  various systems such as superconducting circuits\cite{havlivcek2019supervised} and  photonic systems\cite{cai2015entanglement},
which leads to enormous enthusiasm applying quantum algorithms into various areas\cite{xia2018quantum,hu2020quantum,hu2020quantum,li2021practical,xia2017electronic,sajjan2021quantum,xia2021quantum}.

There is no doubt that we are now in the age of big data and there is an urgent need for developing game-changing quantum algorithms to perform machine learning tasks on large-scale scientific datasets for various industrial and technological applications based on optimization. For a proof of concept,  Du and coworkers have successfully distinguished handwriting numbers '6' and '9' with the quantum support vector machine \cite{li2015experimental}. However, it could be difficult to deal with more challenging problems
, especially when the training data can not be divided apart by a single curve, instead, more than one curve or even enclosed curves might be required to separate them apart.
Another remarkable development  is applying  quantum machine learning   on  variational circuits\cite{schuld2019quantum,arrazola2018machine}, which theoretically, can always be able to classify data with complex distribution. Yet generally, these algorithms rely heavily on a gradient-based systematic optimization of parameters\cite{mitarai2018quantum,farhi2018classification}. On the other hand, quantum nearest neighbor algorithm\cite{wiebe2014quantum} offers  another option to classify data without the gradient based optimization process. In brief, the core of nearest-neighbor classification algorithm is to assign the training vectors into classes, and in each class vectors are close to each other.

	In this study, we will propose a quantum classification algorithm, with which we can build a quantum circuit that is able to classify artificially generated data, and all parameters in the circuit can be obtained without relying on the gradient based optimization process. For this purpose, we introduce 'sublabels' to assist in  classifying  data with intricate distribution, where 'sublabel' represents a minor label subordinates to the main one,  it also called 'subclass'. There are two main tasks in our developed algorithm: how to find the appropriate sublabels and how to build the  quantum classification circuit with these sublabels.
	With the numerical simulation we will demonstrate the application on various classification problems, especially on constructing different phases of materials.
	First, in section 1,  we will present  the basic elements of the algorithm.  Then, we will apply the  algorithm  for classifications for several systems: classification of metallic and insulating phases in the phase diagrame of  VO$_2$;
classification of entanglement in Werner states, and  classification of randomly generated data.
Finally, we will present scaling analysis and discuss generalization  of  quantum classification algorithm in higher dimensional space.
In addition, we present in the supplementary materials all the details of the quantum classification algorithm with examples.

	\section{Algorithm design}
	\label{Algorithm design}

	Consider the training data set $\{{\bf x}_i, y_i\}$, where ${\bf x}_i$ is a vector in $\mathbb{R}^d$, where $d$ is the dimension and $y_i$ represents the label with possible values $\{l_1, l_2, \dots l_{M}\}$. Our goal is to build a quantum classification circuit that can be used to  predict the label for new vectors $\{\bf x_t\}$.
	The classification algorithm is divided into two steps:  The first step is a learning process, where one needs to find the “sublabels” for each class of the training data. Then, based on the information obtained in this learning process, we construct a quantum classification circuit that contains multi control quantum gates.

	In the learning process, firstly we apply the Lloyd's algorithm\cite{mackay2003information} for unsupervised machine learning, which assigns training vectors to the same class as the closest mean vector. However the results derived by Lloyd's algorithm can not be used directly as there might be a sublabel redundancy or not  enough sublabels to reconstruct the initial distribution. To address this issue, in addition to the algorithm for clustering, we propose to use two adjusting algorithms:  one to reduce excessive sublabel and  the other to  make sure that there is no overlap between sublabels.

	For each sublabel, we need to store the  information and build a quantum circuit to estimate the inner product, which is shown in fig.(\ref{device}). When ${\bf x}_i$ are vectors in two-dimensional  space, each data point (or sublable) could be represented by a single qubit.
	An arbitrary state of a single qubit could be written as

	\begin{equation}
	    |\psi\rangle = e^{i\alpha}[\ cos(\theta/2)|0\rangle + sin(\theta/2)e^{-i\phi/2}|1\rangle\ ]
	    \label{1qbstate}
	\end{equation}

	where $\alpha$ is the global phase and the vector ${\bf x}_i= (x_{1,i},x_{2,i})$ is mapped as
	\begin{align*}
	    \theta_i = \frac{2\pi( x_{1,i}-min\{x_1\})}{max\{x_1\}-min\{x_1\}}\\
	    \phi_i = \frac{2\pi (x_{2,i}-min\{x_2\})}{max\{x_2\}-min\{x_2\}}
	\end{align*}
	Here $max\{x_1\}$ and $max\{x_2\}$  represent the maximum value of all $x_{1,i}$  and  $x_{2,i}$ respectively, while $min\{x_1\}$ and $min\{x_2\}$  represent the minimum value of all $x_{1,i}$  and  $x_{2,i}$ respectively.  Then we need to find a measure  describing the distance between the two states, where the 'distance' might be the Euclidean distance between the two vectors\cite{wiebe2014quantum}, or the inner product of their two corresponding quantum states. Here, we chose to calculate the inner product, as calculating the Euclidean distance is more time  and resource-consuming.

	An arbitrary state of a single qubit as shown in Fig.(\ref{1qbstate}),  could be prepared by three rotational  gates:
	\begin{equation}
	    |\psi(\theta_1, \phi_1)\rangle = R_z(\phi_1/2)R_y(\theta_1)R_z(-\phi_1/2)\ |0\rangle
	\end{equation}
	Thus, the inner product between $|\psi(\theta_1, \phi_1)\rangle$ and $|\psi(\theta_2, \phi_2)\rangle$ is given by:
	\begin{equation}
	    \langle\psi(\theta_1, \phi_1)|\psi(\theta_2, \phi_2)\rangle
	    =\langle 0|R_z(\phi_1/2)R_y(-\theta_1)R_z(-\phi_1/2)R_z(\phi_2/2)R_y(\theta_2)R_z(-\phi_2/2)|0\rangle .
	\end{equation}
	The circuit that estimates the inner product will contain six rotational gates, as shown in Fig.(\ref{device}). After a measurement of the final state in the $Z$-basis, the probability of  getting a state  $|0\rangle$ will be an estimation of the inner product $\langle\psi(\theta_1, \phi_1)|\psi(\theta_2, \phi_2)\rangle$. For simplicity, in the following sections, we will write
	\begin{equation*}
	    R(\theta, \phi)=R_z(+\phi/2)R_y(\theta)R_z(-\phi/2)
	\end{equation*}

	\begin{figure}[H]
		\begin{center}
			\centerline{
    \begin{tikzpicture}[thick]
    %

    \tikzstyle{operator} = [draw,fill=white,minimum size=1.5em]
    \tikzstyle{phase} = [fill,shape=circle,minimum size=5pt,inner sep=0pt]
    \tikzstyle{surround} = [fill=blue!10,thick,draw=black,rounded corners=2mm]
    %

    \node at (0,0) (q1) {\ket{0}};
    %

    \node[operator] (op11) at (2,0) {$R_z(\frac{\phi_1}{2})$} edge [-] (q1);
    %

    \node[operator] (op12) at (4,0) {$R_y({-\theta_1})$} edge [-] (op11);
    %

    \node[operator] (op13) at (6,0) {$R_z(-\frac{\phi_1}{2})$} edge [-] (op12);
    %

    \node[operator] (op14) at (8,0) {$R_z(\frac{\phi_2}{2})$} edge [-] (op13);
    %

    \node[operator] (op15) at (10,0) {$R_y({\theta_2})$} edge [-] (op14);
    %

    \node[operator] (op16) at (12,0) {$R_z(-\frac{\phi_2}{2})$} edge [-] (op15);
    \node (end1) at (13.5,0) {} edge [-] (op16);

    \begin{pgfonlayer}{background}
    \node[surround] (background) [fit = (q1) (end1)] {};
    \end{pgfonlayer}
    %
    \end{tikzpicture}
  }
	    \captionsetup{justification=centering}
	    \caption{
				{\bf Sketch of the quantum circuit for estimating the inner product:} Circuit to estimate the  inner product of two-dimensional  vectors contains six rotational  gates. Additionally, we also need a memory to store the information of this group, which could be represented by $\theta_{m}$, $\phi_{m}$, and an integer $N$, which represents the total number of data in this group. }
		\label{device}
		\end{center}
\end{figure}

Moreover,  for every sublabel it is required to store 2 floating  numbers $\theta_m$ and $\phi_m$ that represent the centroid vector of this subgroup, and an integer N representing the total number of data points in this subgroup.

	In the learning process, three  basic algorithms are applied to assist in obtaining “sublabels” from the given training data. Algorithm (\ref{algorithm1}) is designed for an initial  clustering of  the training data. When designing Algorithm (\ref{algorithm1}),  we refer  to Lloyd's algorithm\cite{mackay2003information}in which  we need to assign each vector to the cluster with the closest mean, and then recalculate the centroids of the new cluster. Algorithm (\ref{algorithm1}) will divide the training data with the same prior label into several subgroups.  Algorithm (\ref{algorithm2}) will reduce redundancy, and Algorithm (\ref{algorithm3}) is introduced to make sure there will be no overlap between any two left sublabels of the different prior labels.
    The goal is to leave only the minimal sublabels without losing important information of the training data. After applying Algorithm (\ref{algorithm1}), and repeating Algorithm (\ref{algorithm2}), Algorithm (\ref{algorithm3}) for a number of times, we can get a set with minimal sublabels and information of the centroid vectors for each subgroup (Details of these three algorithms are in the supplementary materials).

	Now, the next step is to build the quantum classification circuit based on the previously obtained information. Consider the following sublabel-control operations,

	\begin{eqnarray*}
	\left.\begin{aligned}
	    &\ {\bf IF} \ SUB LABEL = l_j\\
	    &\ \qquad{\bf DO} \quad ROTATION \ R(-\theta^{l_j}_{m}, -\phi^{l_j}_{m})\\
	    &\ {\bf IF} \ SUB LABEL = l_{j+1}\\
	    &\ \qquad{\bf DO} \quad ROTATION \ R(-\theta^{l_{j+1}}_{m}, -\phi^{l_{j+1}}_{m})\\
	    &\ \cdots\\
	    &\ {\bf FOR}\  SUB LABEL \ in\  PRIOR LABEL \ L_i\\
	    &\qquad{\bf DO} \quad OPERATION \ U(L_i)
	    \end{aligned}\right.
	\end{eqnarray*}

	where the operations $U(L_i)$ are obtained with the  aim to reach the final state $|\Psi^{f}_{L_i}\rangle$. If one wants to set all vectors belong to prior label $L_i$ close to the final state $|\Psi^{f}_{L_i}\rangle$, then $U(L_i)$ can be chosen to satisfy
	\begin{equation*}
	    |\Psi^{f}_{L_i}\rangle = U(L_i)|0\rangle
	\end{equation*}

	To classify a test vector with unknown label, we prefer to rely on a single quantum classification circuit, instead of repeating comparing inner products with the training data. The very basic and intuitive idea is to measure inner product between the new one and all centroid vectors of subgroups. Circuits for classification will consist of two parts: the control qubits representing the sublabels, and others representing the given new vector. Fig.(\ref{device2}) is a sketch showing the structure of the main circuit. First, one map the test data ${\bf x_t}$ into the prepared circuit as
	\begin{equation}
	    {\bf x_t}\rightarrow|\Psi({\bf x_t})\rangle=\frac{1}{\sqrt{n}}\left[\sum_i^n|l_i\rangle\right]\bigotimes|\psi(\theta,\phi)\rangle
	\end{equation}
	where

	$|l_i\rangle$ are the orthogonal eigenstates for the control qubits, representing sub labels $l_i$, and there are $n$ sublabels in total. In the mapping process, we need to apply Hadamard gates for L-qubits (representing the sublabels), and apply operator $T({\bf x_t})$ on the V-qubits (representing the test data), where $|\psi(\theta,\phi)\rangle=T({\bf x_t})|0\rangle$. The quantum classification circuit can be described as:
	\begin{equation}
	    U = \sum_i^n\left[|l_i\rangle\langle l_i|\bigotimes U^{}(L_i) R(-\theta^{l_i}_{m}, -\phi^{l_i}_{m})\right] + \sum_{j=n+1}^{2^N-1}\left[|j\rangle\langle j|\bigotimes I\right]
	\end{equation}
	where there are $N$ qubits used to represent sublabels, and $n$ sublabels totally. We can notice that for $0\leq k\leq n$,
	\begin{equation}
	    \begin{split}
	        \langle l_k, \Psi^{f}_{L_k}| U|\Psi({\bf x})\rangle
	        &=\langle l_k, \Psi^{f}_{L_k}|\sum_i^n
	        \left[|l_i\rangle\bigotimes U^{}(L_i) R(-\theta^{l_i}_{m}, -\phi^{l_i}_{m})|\psi(\theta,\phi)\rangle\right]
	        +\langle l_i, \Psi^{f}_{L_k}|\sum_{j=n+1}^{2^N-1}\left[|j\rangle\langle j|\bigotimes I\right]
	        \\
	        &=\sum_i^n\left[
	        \delta_{ik} \langle\Psi^{f}_{L_i}|U^{}(L_i) R(-\theta^{l_i}_{m}, -\phi^{l_i}_{m})|\psi(\theta,\phi)\rangle
	        \right]
	        \\
	        &=
	        \langle\psi(-\theta^{l_k}_{m}, -\phi^{l_k}_{m})|\psi(\theta,\phi)\rangle
	    \end{split}
	\end{equation}
	By applying a measurement with eigenstate $|l_k\rangle|\Psi^{f}_{L_k}\rangle$, the inner product $\langle\psi(-\theta^{l_k}_{m}, -\phi^{l_k}_{m})|\psi(\theta,\phi)\rangle$ can be estimated.

	Fig.(\ref{device2}) is a sketch of the main circuit, where $T({\bf x_t})$ represents the given test data, $l$ is the sublabel and $L$ is prior label.

	\begin{figure}[H]
			\centerline{
    \begin{tikzpicture}[thick]
    \tikzstyle{operator} = [draw,fill=white,minimum size=1.5em]
    \tikzstyle{phase} = [fill,shape=circle,minimum size=5pt,inner sep=0pt]
    \tikzstyle{surround} = [fill=blue!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_blue} = [fill=blue!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_red} = [fill=red!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_yellow} = [fill=yellow!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_green} = [fill=green!10,thick,draw=black,rounded corners=2mm]
    %

    \node at (0,0) (q1) {L-qubits\qquad\ket{0}};
    \node at (0,-1) (q2) {V-qubits\qquad\ket{0}};
    %
    \node[operator] (op11) at (2,0) {$H^{\otimes N}$} edge [-] (q1);
    \node[operator] (op21) at (2,-1) {$T(\bf x_t)$} edge [-] (q2);

    \node[circle0] (op12) at (4.5,0) {$l_1$} edge [-] (op11);
    \node[operator] (op22) at (4.5,-1) {$U(L_{l_1})R(-\theta_1,{-\phi_1})$} edge [-] (op21);
    \draw[-] (op12) -- (op22);
    \node[] (achor0) at (3,0.25) {};
    \node[circle0] (op13) at (8,0) {$l_2$} edge [-] (op12);
    \node[operator] (op23) at (8,-1) {{$U(L_{l_2})R(-\theta_2,{-\phi_2})$}} edge [-] (op22);
    \draw[-] (op13) -- (op23);
    \node[circle0] (op14) at (11.5,0) {$l_n$} edge [-] (op13);
    \node[operator] (op24) at (11.5,-1) {{$U(L_{l_n})R(-\theta_n,{-\phi_n})$}} edge [-] (op23);
    \draw[-] (op14) -- (op24);
    \node[meter] (meter1) at (14,0) {} edge [-] (op14);
    \node[meter] (meter2) at (14,-1) {} edge [-] (op24);
    \draw [dashed] (9.75,0.75) -- (9.75,-1.75);
    \begin{pgfonlayer}{background}

    \node[sd_green] (background) [fit = (op11) (op21)] {};
    \node[sd_blue] (background) [fit = (achor0) (op24)] {};
    \node[sd_red] (background) [fit = (meter1) (meter2)] {};
    \end{pgfonlayer}
    %
    \end{tikzpicture}
  }
	    \captionsetup{justification=centering}
	    \caption{
				{\bf Sketch of the main circuit :} Qubits in the main circuit can be divided into two groups: L-qubits to represent the sub labels, and will play the role of control bits and the  V-qubits to represent a  given vector. Initially, L-qubits will be prepared at state $(\frac{|0\rangle+|1\rangle}{\sqrt{2}})^{\otimes N}$, where there are N qubits representing the sub labels. The minimum $N = \lceil\log_2 n\rceil$, and n is the total number of sub labels. In sum $n$ control rotation gates are needed.}
		\label{device2}
	\end{figure}

	When predicting the label for a new test vector, we assume that the probability for each possible sublabels is the same. Based on this assumption, we applied $N$ Hadamard gates on the label qubits, preparing them as a uniform superposition state. In fact, one can always adjust the probability for sublabels based on the training result, and only keep the states representing 'valid' sublabels. Assume that we finally derived $K$ sublabels, and $N$ qubits are used as label qubits, where $2^{N-1}<K\leq 2^N$.
	For a training data set with two labels $\{{\bf x_i}, y_i\}$, $y_i\in\{0, 1\}$, $K_0$ sublabels are assigned with label $y=0$, while the other $K_1$ labels with $y=1$, $K_0+K_1=K$.
	Then the Hadamard gates $H^{\otimes N}$ can be replaced as some certain operation, in order to prepare the label qubits at state
	\begin{equation}
	    |\phi_L\rangle=
	    \sum_{n=0}^{K_0-1}\sqrt{p_n}|n\rangle +
	    \sum_{m=0}^{K_1-1}\sqrt{p_{m+K_0}}|m+2^{N-1}\rangle
	    \label{lqubits}
	\end{equation}
	where $p_n\geq0$, and $\sum_{n=0}^{K-1}p_n=1$. Obviously, for sublabels assigned with label $y=0$, the first label qubit is set as state $|0\rangle$, and for the others the first qubit is set as $|1\rangle$. For a new test data, one only needs to measure the first label qubit and the data qubits. By comparing $P(q_1=|y\rangle, V=|\Psi^{f}_{y}\rangle)$, $y_i\in\{0, 1\}$, the new test vector will be assigned with label corresponding to the larger $P$, where $P(q_1=|y\rangle, V=|\Psi^{f}_{y}\rangle)$ represents probability to find the first label qubit $q_1$ at state $|y\rangle$ and meanwhile find the data qubits at state $|\Psi^{f}_{y}\rangle$. $p_n$ are chosen to maximize
	\begin{equation}
	    \sum_i\left\{
	    P(q_1=|y_i\rangle, V=|\Psi^{f}_{y_i}\rangle)
	    -\lambda P(q_1\neq|y_i\rangle, V=|\Psi^{f}_{y_i}\rangle)
	    \right\}
	\end{equation}
	where $\lambda\geq0$ is the penalty.

	\section{Applications}
	\label{Applications}
	\subsection{Classification of metallic-insulating  phases of  vanadium dioxide}
	\label{Classification of metallic-insulating  phases of  vanadium dioxide}
	Strongly correlated electron materials and their phase transitions have attracted great interest for device application and  in condensed matter physics\cite{dagotto2005complexity,arko1997strongly}. Recently, metallic-insulating phase transition of vanadium dioxide ($VO_2$), as a prototype of strongly correlated electronic materials, attracted experimental and theoretical attention for its distinct structures and electronic properties \cite{qazilbash2007mott,jeong2013suppression}.

Here, we will apply the developed quantum classification algorithm to distinguish the metallic state from the insulating  state of $VO_2$. Data used in this section is based on experimental results reported in Ref. \cite{chen2017pressure}.  $VO_2$ exhibits several special structures under different temperature and pressure. As shown in Fig.(\ref{subfig_a_vo2}), the red dots represent metallic state, the blue dots represent the insulating state, and the black solid line represents the phase transition line. Note, our training data were chosen far from the phase transition line in order to test the classification power of the designed quantum algorithm.  Fig.(\ref{subfig_b_vo2}) show the initial clustering results after applying the Algorithm (\ref{algorithm1}) once. 	In Algorithm (\ref{algorithm1}) one should determine the parameter $D\in(0,1]$, which is the minimum inner product between the centroid vector and arbitrary vectors with the same sublabel. In other words, a vector can be assigned to a certain sublabel when the inner product between itself and the centroid vector of the sublabel are larger than $D$.  Here we set $D=0.99$. Now, we  need to repeat Algorithm (\ref{algorithm2}) and Algorithm (\ref{algorithm3}) several times to reduce the excessive sublabels as shown in Fig.(\ref{subfig_c_vo2}).  After repeating Algorithm(\ref{algorithm2}) and Algorithm (\ref{algorithm3}) three times, the number of classes can not be reduced further as shown
in Fig.(\ref{subfig_c_vo2}).

In both Fig.(\ref{subfig_b_vo2}) and Fig.(\ref{subfig_c_vo2}),  blue and red spheres are used to represent data with the same sublabel, where the center of sphere represents the average vector, and their radius is used to represent the number of vectors belonging  to a  certain sublabel. In Fig.(\ref{subfig_d_vo2}), we demonstrate the prediction results of arbitrary given data. States in the yellow parts are predicted to behave  as metallic states, and the  ones in the blue parts are predicted to be  insulating states. In the training process, we used 100 data vectors for the metallic states and 100 data vectors for the insulating states, out of the totally 1100 data vectors.

	Simulation results shows that our developed quantum algorithm can efficiently classify metallic or insulating  states of $VO_2$.   Finally we need 7 sublabels for insulating  states and 8 sublabels for metallic states, which means that the classification circuit consists of 5 qubits (4 for sublabels as control qubit and 1 for data).

It is important to note,  when Pressure ($P$) or  Temperature ($T$) is small, prediction can be incorrect. The error appears because of the fact that few vectors in these area where  used in the training process. Moreover, when we convert a vector $(P, T)$ into quantum states, we changed them as $\theta, \phi\in[0,2\pi)$. Though classically $T= 0{}^oC$ and $T= 120{}^oC$ are extremely different, when we convert them as angles, $\theta=0$ and $\theta=2\pi$ are nearly the same. One might notice that there is a slim yellow line around $P=20GPa$ in our prediction. It is reported in Ref. \cite{chen2017pressure} that there is a structural transition  between the state $M1$ and $M1'$. However, this is not a metallic-insulating  transition. In our classification algorithm, we did not expect to predict this transition and the slim line shows up “coincidentally”. To get a better prediction results, an option is to map both training data and test data into $[0, \pi]\times[0, \pi]$ instead of $[0, 2\pi]\times[0, 2\pi]$. However, here we do not focus on the  performance at low temperature, nor the phase transition about $P=20GPa$. As a result, mapping data into $[0, 2\pi]\times[0, 2\pi]$ is still an acceptable choice. More simulation results of the phase transition with different training data are offered in the supplementary materials.

	\begin{figure}[H]

\centering

\subfigure[]{
\centering
\label{subfig_a_vo2}
\includegraphics[width=0.45\linewidth]{fig_data_vo2.pdf}
}
\subfigure[]{
\centering
\label{subfig_b_vo2}
\includegraphics[width=0.45\linewidth]{fig_sublabel_vo2.pdf}
}

\vfill

\subfigure[]{
\label{subfig_c_vo2}
\includegraphics[width=0.45\linewidth]{fig_reduce_vo2.pdf}
}
\hspace{0.01\linewidth}
\subfigure[]{
\label{subfig_d_vo2}
\includegraphics[width=0.45\linewidth]{classify_vo2.pdf}
}
\caption{{\bf Classification of  metallic  and insulating  states of $VO_2$}: (a)Initial data used for classification. Red dots represent metallic state, and blue ones represent insulating state. Phase transition line indicated by  the black solid curve.(b) Forming subgroups after applying algorithm[\ref{algorithm1}] once. Similarly in (b) and (c), blue or red spheres are used to represent data with the same sublabel, where the center of sphere represents the average vector, and the radius represents number of vectors belong to this sublabel. (c) Results after repeating algorithm[\ref{algorithm2},\ref{algorithm3}] 3 times. (d)Prediction of new data. New vector in the blue part will be recognized with label 'insulating', and label of new vectors in yellow part will be predicted as 'metallic'. Blue and red dots are still the initial data.}
\label{fig_simu_vo2}
\end{figure}

    \subsection{Classification of randomly  generated data}
    \label{Classification of randomly  generated data}

Here we will show another classification example, where the distribution of training data is artificially generated randomly.
Different from the example of $VO_2$, here the two groups can not be divided with a simple single boundary.
We generated 1100 red points and 1100 blue points at random, from which 100 red points and 100 blue points are picked up randomly as training data, the others will be used as test data.
    All test points are shown in Fig.(\ref{subfig_a_complex}), and the training points in Fig.(\ref{subfig_b_complex}).
    The two isolated groups of red points in  Fig.(\ref{subfig_a_complex}) are scattered along with the  blue ones covering  the whole area as shown in Fig.(\ref{subfig_b_complex}). The distribution of initial data makes it more challenging to classify the test data. After appropriate learning process, 54 sublabels (22 for red and 32 for blue) are obtained from the training data. Finally a classification circuit can be build with 7 qubits, 6 as control qubits and 1 for the given data. Prediction for labels of test data is shown in Fig.(\ref{subfig_c_complex}), where light blue dots are training points that are predicted as 'blue', yellow dots are predicted as 'red', red and blue cross represent the training data. Totally, 878 red test data and 836 blue test data are classified correctly, the matching rate for red and blue points can be estimated  as $87.8\%$ and $83.6\%$ respectively (when calculating match rate the training data are all excluded).

    \begin{figure}[H]

\centering

\subfigure[]{
\centering
\label{subfig_a_complex}
\includegraphics[width=0.31\linewidth]{complex_test.pdf}
}
\subfigure[]{
\centering
\label{subfig_b_complex}
\includegraphics[width=0.31\linewidth]{complex_training.pdf}
}
\subfigure[]{
\label{subfig_c_complex}
\includegraphics[width=0.31\linewidth]{complex_classify.pdf}
}

\caption{{\bf Classification of randomly generated data:}
(a) Training data. In total, we generated 1100 blue points and 1100 red pints by the same generating function. 100 red points and 100 blue points are chosen randomly as training data, and the left pints are used as test data. (b) Training data including 100 red points are 100 blue points. (c) Prediction for the labels of test data. Light blue points are test data predicted as 'blue', and yellow points are predicted as 'red'. Red and blue cross represent training data.}
\label{fig_simu_complex}
\end{figure}

    \subsection{ Entanglement classification  in  Werner states}
    \label{Entanglement classification  in  Werner states}
    Further, our method can also be applied in entanglement classification.
    Consider the following scenario, Alice wants to send a  message to Bob, in which she will send  some entangled photon pairs of states  as digit 1 and some pairs at an untangled state as digit 0. Initially, Alice will send some photon pairs at various states for training by  informing   Bob which pair represents 1 and which represents 0.  Later she will use photon pairs states for communication. Although the widely used CHSH inequality\cite{clauser1969proposed} can be used to detect entanglement as the violation of CHSH inequality guarantees the existence of entanglement, however we can’t make any conclusion if the inequality is not violated. To address this issue consider Werner states in the density matrix form:
    \begin{equation}
        \rho_W(p,\phi)=\frac{p}{4}|\Psi_B(\phi)\rangle\langle\Psi_B(\phi)|+\frac{1-p}{4}I
        \label{werner}
    \end{equation}
    where $I$ is  the $4\times4$  identity matrix, the parameter $0<p<1$, and   $|\Psi_B(\phi)\rangle$ is the Bell state given by:
    \begin{equation}
        |\Psi_B(\phi)\rangle=\frac{1}{\sqrt{2}}\left(
        |\uparrow\downarrow\rangle+e^{i\phi}|\downarrow\uparrow\rangle
        \right)
    \end{equation}

    Assume that Alice uses Werner states to transport information, while Bob will carry on the Bell test experiment with the following measurements: $Z, X; \frac{Z+X}{\sqrt{2}}$, and $ \frac{Z+X}{\sqrt{2}}$. From the measurement results, Bob will calculate four important correlation functions $E(Z, \frac{Z+X}{\sqrt{2}}), E(X, \frac{Z+X}{\sqrt{2}}), E(Z, \frac{Z-X}{\sqrt{2}}), E(X, \frac{Z-X}{\sqrt{2}})$ where
    \begin{equation}
        E(a,b)=\frac{N_{++}+N_{--}-N_{+-}-N_{-+}}{N_{++}+N_{--}+N_{+-}+N_{-+}}
    \end{equation}
    $N_{++}$ are the number of photon pairs whose measurement results are both $+1$ in the two channels. If Alice sets $\phi=0,\pi$, Bob will observe violation of CHSH inequality for $p>\frac{1}{\sqrt{2}}$. If Alice sets  $\phi=\pm\frac{\pi}{2}$, Bob can never observe violation of CHSH inequality. However, $\rho_W(p,\phi)$ will be entangled state when $p>\frac{1}{3}$. Consequently, CHSH inequality will not be a good classification way  for Bob. Instead, if Bob can set up a machine learning based on neural network, he will be able to 'decode' Alice's information with a much higher match rate\cite{gao2018experimental}. Here, we will show that our quantum classification algorithm can classify entanglement states  in such Werner state. We will take the 4-dimesinal vectors $E(Z, \frac{Z+X}{\sqrt{2}}), E(X, \frac{Z+X}{\sqrt{2}}), E(Z, \frac{Z-X}{\sqrt{2}})$, and  $E(X, \frac{Z-X}{\sqrt{2}})$ as input calculated  based on measurement results.

    By changing $\phi$ and $p$ in Eq.(\ref{werner}), we can generate different entangled or untangled states. We prepared 200 entangled states and 200 untangled states as the test data. Moreover, we also generated a few different training data set, in each set half are entangled and the other half are untangled. After learning based on different training set, we can build quantum classification circuit to distinguish entanglement from the test data, and the simulation results are shown in Fig.(\ref{fig_simu_werner}). In (a), the training set only contains 32 points, and we keep them all as sublabels. So that there are $7$ qubits in the classification circuit (5 for sublabels and 2 for test data). 12 points are predicted with wrong label. In (b) the training set contains 64 points and all are kept as sublabels, and there are 8 qubits in the classification circuit. Finally 8 points are predicted with wrong labels. In (c) there are 128 points in the learning set and we derived 64 sublabels. Similarly as (b) we need 8 qubits to build the classification circuit, and only 6 qubits are predicted with wrong sublabels. However in (c) we do not plot the sublabels, as by applying our classification algorithms we can only get the parameters $E$ for every sublabel instead of $r, \phi$.

    \begin{figure}[H]

\centering

\subfigure[]{
\centering
\label{subfig_a_werner}
\includegraphics[width=0.45\linewidth]{werner_classify16_8.pdf}
}
\subfigure[]{
\centering
\label{subfig_b_werner}
\includegraphics[width=0.45\linewidth]{werner_classify32_6.pdf}
}
\subfigure[]{
\label{subfig_c_werner}
\includegraphics[width=0.45\linewidth]{werner_classify64_3.pdf}
}

\caption{{\bf Entanglement classification for Werner states:} In the plots, the radius represents the parameter $p$, and the angle represents  $\phi$. Every single dot represents a Werner state. Yellow dots represent test data that are predicted as 'untangled', and light blue dots represent test data predicted as 'entangled'. Cross represent training data, and red cross for untangled states, blue for entangled states. In all three figures we used the same 400 test data, half of which are entangled and the other half are untangled. Half of the training data are entangled states and the other half are untangled states. (a)32 test data are used, all are kept as sublabels. (b)64 test data are used, all are kept as sublabels. (c)128 test data are used, and only 64 of them are used as sublabels.}
\label{fig_simu_werner}
\end{figure}

    In these figures the points $(r,\phi)$ represent  Werner states.  Notice that $r, \phi$ are not used in the learning or classification process, as Bob does not know exact $r, \phi$ either. In the supplementary materials we  provide details of  the simulation.

    In the above discussion we assume that in communication between Alice and Bob, all measurement results are discrete, and one can easily calculate the parameters $E$. However, in chemical reactions the measurement results are often continuous, and one will get some special distributions after measurement. Recently, Zare and coworkers \cite{perreault2018cold} reported the rotationally inelastic collisions of deuterium hydride (HD, prepared at certain quantum states) with H$_2$ and D$_2$ under extremely low temperature (mean collision energy around 1 K), and they found that the orientation of HD molecules will lead to different distribution of scattering angle\cite{perreault2018cold}.
    If the scattering experiment is applied as measurement to detect entanglement, then it would be nearly impossible to derive information about entanglement directly from the intricate raw data. Under such situations, some special methods, as we discussed in ref.\cite{li2019entanglement}, are required. With the assistance of auxiliary functions it will be possible to obtain the parameters $E$ from raw measurement results, after which by the same procedure we can build a quantum circuit for classification.

	\section{Further Discussions}
	\label{Further Discussions}
	 So far, our classification has been restricted to two and four-dimensional vectors.Here, we discuss how to use it to classify vectors in higher dimensional space. Depending on the structure of qubits, two different mapping methods can be used: {\bf Mapping method I:} An arbitrary quantum state of N-qubits can be described as
	\begin{equation*}
	    |\Psi\rangle= \sum_{i=0}^{2^N-1}{c_i|i\rangle}
	\end{equation*}
	where $c_i$ can be rewritten as a function of ${\bf \Theta}=(\theta_1,\cdots,\theta_{2^N-1},\phi_1,\cdots,\phi_{2^N-1})$:
	\begin{align*}
	    &c_0 = \cos\theta_1\\
	    &c_1 = e^{i\phi_1}\sin{\theta_1}\cos{\theta_2}\\
	    &\cdots\\
	    &c_{2^N-2} = e^{i\phi_{2^N-2}}\Pi_{j=1}^{2^N-2}\sin{\theta_j}\cos{\theta_{2^N-2}}\\
	    &c_{2^N-1} = e^{i\phi_{2^N-1}}\Pi_{j=1}^{2^N-1}\sin{\theta_j}
	\end{align*}
	and $0\leq\theta_j\leq\pi$,$0\leq\phi_j\leq 2\pi$. Then there exists a  mapping operator $T(\bf \Theta)$, $|\Psi({\bf \Theta})\rangle=T(\bf \Theta)|0\rangle$, by which one could map a vector in $(2^N-1)$-dimensional  space into a quantum state of N-qubits. For this mapping method, the main structure of circuit is still like the one shown in Fig.(\ref{device2}). If all qubits in the main circuit are connected with others and we could build arbitrary quantum gates between any qubits in the main circuit, then we can obviously apply this mapping method I. However, sometimes connection in the machine could not satisfy our demand, then the second  mapping method might be more acceptable {\bf Mapping method II:} An arbitrary untangled quantum state of N-qubits can be described as:
	\begin{equation*}
	    \bigotimes_{j=1}^{N}{T_j({\bf x_t})|0\rangle} = \bigotimes_{j=1}^{N}{\left[\cos{\theta_j|0\rangle}+e^{i\phi_j}\sin{\theta_j}|1\rangle\right]}
	\end{equation*}
	Then we could map the vector ${\bf \Theta}=(\theta_1,\cdots,\theta_{N},\phi_1,\cdots,\phi_{N})$ into the untangled quantum state. Method II requires a circuit where the qubits representing sub labels are connected with all the qubits representing our vector, yet the connection between the 'data qubits' are not required. A sketch of the main circuit using method II can be found in fig.(\ref{device4}), where for simplicity, we note that $U_{ij}=U_j^{}(L_i) R_j(-\theta^{l_i}_{m}, -\phi^{l_i}_{m})$

	\begin{figure}[H]
			\centerline{
    \begin{tikzpicture}[thick]
    \tikzstyle{operator} = [draw,fill=white,minimum size=1.5em]
    \tikzstyle{phase} = [fill,shape=circle,minimum size=5pt,inner sep=0pt]
    \tikzstyle{surround} = [fill=blue!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_blue} = [fill=blue!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_red} = [fill=red!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_yellow} = [fill=yellow!10,thick,draw=black,rounded corners=2mm]
    \tikzstyle{sd_green} = [fill=green!10,thick,draw=black,rounded corners=2mm]
    %

    \node at (0,0) (q1){L\ \ket{0}};
    \node at (0,-1) (q2){V1\ \ket{0}};
    \node at (0,-2) (q3){V2\ \ket{0}};
    %
    \node[operator] (op11) at (1.5,0) {$H^{\otimes N}$} edge [-] (q1);
    \node[operator] (op21) at (1.5,-1) {$T_1(\bf x_t)$} edge [-] (q2);
    \node[operator] (op31) at (1.5,-2) {$T_2(\bf x_t)$} edge [-] (q3);

    \node[circle0] (op12) at (2.75,0) {$l_1$} edge [-] (op11);
    \node[operator] (op22) at (2.75,-1) {$U_{11}$} edge [-] (op21);
    \draw[-] (op12) -- (op22);
    \node[circle0] (op13) at (4,0) {$l_1$} edge [-] (op12);
    \node[operator] (op33) at (4,-2) {$U_{12}$} edge [-] (op31);
    \draw[-] (op13) -- (op33);
    \node[circle0] (op14) at (5.25,0) {$l_2$} edge [-] (op13);
    \node[operator] (op24) at (5.25,-1) {$U_{21}$} edge [-] (op22);
    \draw[-] (op14) -- (op24);
    \node[circle0] (op15) at (6.5,0) {$l_2$} edge [-] (op14);
    \node[operator] (op35) at (6.5,-2) {$U_{22}$} edge [-] (op33);
    \draw[-] (op15) -- (op35);
    \node[circle0] (op16) at (7.75,0) {$l_3$} edge [-] (op15);
    \node[operator] (op26) at (7.75,-1) {$U_{31}$} edge [-] (op24);
    \draw[-] (op16) -- (op26);
    \node[circle0] (op17) at (9,0) {$l_3$} edge [-] (op16);
    \node[operator] (op37) at (9,-2) {$U_{32}$} edge [-] (op35);
    \draw[-] (op17) -- (op37);
    \node[meter] (meter1) at (10.5,0) {} edge [-] (op17);
    \node[meter] (meter2) at (10.5,-1) {} edge [-] (op26);
    \node[meter] (meter3) at (10.5,-2) {} edge [-] (op37);
    \begin{pgfonlayer}{background}
    \node[sd_green] (background1) [fit = (q1) (op31)] {};
    \node[sd_blue] (background2) [fit = (op12) (op37)] {};
    \node[sd_red] (background3) [fit = (meter1) (meter3)] {};
    \end{pgfonlayer}
    %
    \end{tikzpicture}
  }
	    \captionsetup{justification=centering}
	    \caption{
				{\bf Sketch of the main circuit for Method II:}
				Qubits in the main circuit can be divided into two groups: One group will represent the sub labels, and will play the role of control bits, as the $L$ part in the figure. The other will represent the given vector, as the $V$ part in the circuit. Furthermore, qubits representing the vectors are divided into a few groups (In this Figure, 2 groups), and the sublabel qubits will control them respectively. We need to measure all of them to get our results. Connection between the $V$ qubits are not required in this circuit.}
		\label{device4}

	\end{figure}

	For the complexity analysis, let us assume that M d-dimensional vectors are offered as training data, thus $\lceil \log_2{d}\rceil$ qubits are required to represent the vectors. When measuring the inner product of two vectors, $\mathcal{O}\left\{\exp{(\log_2{d})}\right\}$ times of measurements are required. In the learning process to obtain sublabels, as discussed in sec.I, we need to repeat calculating inner products for $\mathcal{O}\left (M^2\right)$ times. Totally, the time complexity to obtain sublabels is $\mathcal{O}\left (M^2d\right)$. Then assume that we finally obtained  $L$ sublabels, then $\lceil \log_2{L}\rceil$ qubits are required to represent the sublabel. Thus, the quantum circuit to predict labels of test data contains $L$ multi-control gates. If we prepare the label qubits at uniform superposition state, then all qubits should be measured at last. For the label qubits, $\mathcal{O}\left\{\exp{(\log_2{L})}\right\}$ times of measurements are required, while for the data qubits $\mathcal{O}\left\{\exp{(\log_2{d})}\right\}$ times of measurements are required. As a result, we need to repeat measurement for $\mathcal{O}\left (Ld\right)$ times. However, if we prepare the label qubits at state eq.(\ref{lqubits}), then we only need to measure the first label qubit, and required times for measurement will be $\mathcal{O}\left (d\right)$.

	\section{Conclusion}
	\label{Conclusion}

	In summary, we developed a quantum classification algorithm where the training data is firstly clustered and assigned as various sublabels, and then based on these sublabels the quantum circuit is built for classification.
	Further we applied this method for classifications of metallic-insulating  transition in $VO_2$, distinguish entanglement in Werner states, and classify some randomly generated data.
	Numerical simulation result shows that our algorithm is capable for various classification problems, especially the study of phases transition in materials.

\section*{Acknowledgement}
The authors would like to thank Dr. Masoud Mohseni, Dr. Manas Sajjan and Dr. Zixuan Hu for the helpful suggestions and discussions.
We acknowledge the financial support in part by the National Science Foundation under award number 1955907 and funding by the U.S. Department of Energy (Office of Basic Energy Sciences) under Award No.de-sc0019215.