\section{Introduction}\label{sec:introduction}
Deep learning methods can achieve outstanding performance in multiple fields including computer vision~\cite{KXSJ16}, natural language processing~\cite{LFMK15,DKY14}, speech and audio processing~\cite{MN17}, and graph analysis~\cite{JiaweiZhang2019Tutorial}. Training deep learning models involves an optimization process to find the parameters in the model that minimize the loss function. At the same time, the number of parameters commonly used in deep learning methods can be very huge.

Therefore, optimization algorithms are critical for deep learning methods: not only the model performance, but also training efficiency are greatly affected.
In order to cope with the high computational complexity of training deep learning methods, stochastic gradient descent (SGD)~\cite{sgd_overview} is utilized to update parameters based on the gradient of each training sample instead. The idea of momentum~\cite{moment}, inspired by Newton's first law of motion, is used to handle the oscillations of SGD. SGD with momentum~\cite{IJGG13} achieves the faster convergence rate and better optimization results compared with the original SGD. In gradient descent based optimization, training efficiency is also greatly affected by the learning rate.
AdaGrad~\cite{adagrad} is the first optimization algorithm with adaptive learning rates, which makes use of the learning rate decay. AdaDelta~\cite{adadelta} subsequently improves AdaGrad to avoid the extremely small learning rates. {\textsc{ADAM}}~\cite{ADAM} involves both adaptive learning~\cite{sgd_overview} and momentum~\cite{moment} and utilizes the exponential decay rate $\beta_1$ (momentum weight) to accelerate the convergence in the relevant directions and dampen oscillations. However, the decay rate $\beta_1$ of the first-order momentum $\mathbf{m}_t$ in {\textsc{ADAM}} is a fixed number, and the selection of the hyperparameter $\beta_1$ may affect the performance of {\textsc{ADAM}} greatly. Commonly, $\beta_1 = 0.9$ is the most widely used parameter as introduced in~\cite{ADAM}, but there is still no theoretical evidence proving its advantages.

During the optimization process, it is common that there exist errors in some update steps. These errors can be caused by the inappropriate momentum calculation, and then lead to slower convergence or oscillations. For each parameter updates, using the fixed momentum weight fails to take the different influence of the current gradient into consideration, which will render errors in momentum computing. For example, when there exist parts of opposite eigen components~\cite{moment} between the continuous two parameter updates (we regard this situation as an error), the current gradient should be assigned a larger weight to correct the momentum in the last update, instead of being placed with a fixed influence.
We will illustrate this problem through cases in Section~\ref{subsubsec: beta_motivation} where {\textsc{ADAM}} with a fixed weight $\beta_1$ cannot handle some simple but intuitive convex optimization problems. Based on this situation, we need to control the influence of momentum by an adaptive weight. What's more, designing hyperparameter-free optimization algorithms has been a very important research problem in recent years, controlling the number of hyperparameters will not only stabilize the performance of the optimization algorithm, but also release the workload of hyperparameters tuning.

In this paper, we introduce a novel optimization algorithm, namely {\text{DEAM}} (\textbf{D}iscriminative w\textbf{E}ight on \textbf{A}daptive \textbf{M}omentum) to deal with the aforementioned problems. {\text{DEAM}} proposes an adaptive momentum weight $\beta_{1,t}$, which will be learned and updated in each training iteration automatically.
Besides, {\text{DEAM}} employs a novel backtrack term $d_t$, which will restrict redundant updates when {\text{DEAM}} decides that the correction of the previous step is needed.
We also provide the theoretic analyses about the adaptive momentum weight along with extensive experiments. Based on them, we verify that the adaptive momentum term weight $\beta_{1,t}$ and the operation of backtrack term $d_t$ can be crucial for the performance of the learning algorithms.

Here, we summarize the detailed learning mechanism of {\text{DEAM}} as follows:
\begin{itemize}
	\item {\text{DEAM}} computes adaptive momentum weight $\beta_{1,t}$ based on the ``discriminative angle" $\theta$ between the historical momentum and the newly calculated gradient.
	\item {\text{DEAM}} introduces a novel backtrack term, i.e., $d_t$, which is proposed to correct the redundant update of the previous training epoch if it is necessary.
	The calculation of $d_t$ is also based on the discriminative angle $\theta$.
	\item {\text{DEAM}} involves fewer hyperparameters than the {\textsc{ADAM}} during the training process, which can decrease the workload of hyperparameter tuning.
\end{itemize}

Detailed information about the learning mechanism and the concepts mentioned above will be described in the following sections. This paper will be organized as follows. In Section 2, we will cover some related works about widely used optimization algorithms. In Section 3, we will analyze more detail of our proposed algorithm, whose theoretic convergence rate will also be studied. Extensive experiments will be exhibited in Section 4. Finally, we will give a conclusion of this paper in Section 5.
\section{Related Works} \label{sec:relatedwork}
\noindent\textbf{Stochastic Gradient Descent}: Stochastic gradient descent (SGD)~\cite{sgd_overview,bianchi2012convergence} performs parameter updating for each training example and label. The advantages of SGD include fast converging speed compared with gradient descent and preventing redundancy~\cite{sgd_overview}.~\cite{on_variance_reduction} use the variance reduction methods to accelerate the training process of SGD.

\noindent\textbf{Adaptive Learning Rates}: To overcome the problems brought by the unified learning rate, some variant algorithms applying adaptive learning rate~\cite{behera2006adaptive} have been proposed, such as AdaGrad~\cite{adagrad}, {\text{AdaDelta}}~\cite{adadelta}, RMSProp~\cite{rmsprop}, {\textsc{ADAM}}~\cite{ADAM} and recent ESGD~\cite{dauphin2015equilibrated}, AdaBound~\cite{luo2019adaptive}. {\text{AdaGrad}} adopts different learning rates to different variables. One drawback of AdaGrad is that with the increasing of iteration number $t$, the adaptive term may inflate continuously, which leads to a very slow convergence rate in the later stage of the training process. {\text{RMSProp}} can solve this problem by using the moving average of historical gradients.

\noindent\textbf{Momentum}: Momentum~\cite{moment,impor_mom_deeplearning,li2017convergence,dozat2016incorporating,mitliagkas2016asynchrony} is a method that helps accelerate SGD in the relevant direction and prevent oscillations on the descent route. The momentum accelerates updates for dimensions whose gradients are in the same direction as historical gradients, and decelerates updates for dimensions whose gradients are the opposite. Momentum is also applied to Nesterov accelerated gradient (NAG)~\cite{nag}, which provides the momentum term with the estimating next position capability. Instead of using the current location to calculate the gradient, NAG first approximates the next position of the variables, then uses the approximated future position to compute the gradient.
{\textsc{ADAM}}~\cite{ADAM,normalized_direction} is proposed based on momentum and adaptive learning rates for different variables. {\textsc{ADAM}} records the first-order momentum and the second-order momentum of the gradients using the moving average, and further computes the bias-corrected version of them. Based on {\textsc{ADAM}}, Keskar et al.~\cite{switch_from_adam_sgd} proposes to switch from {\textsc{ADAM}} to SGD during the training process. {\text{AMSGrad}}~\cite{amsgrad} is a modified version of {\textsc{ADAM}}, which redefines second-order momentum by a maximum function.

\section{Proposed Algorithm}\label{sec:method}
Our proposed algorithm {\text{DEAM}} is presented in Algorithm~\ref{alg:cwm}. In the algorithm, $f_1, f_2,\dots, f_T$ is a sequence of loss functions computed with the training mini-batches in different iterations (or epochs).
\begin{algorithm}[tb]

	\footnotesize
	\caption{{\text{DEAM}} Algorithm}
	\hspace*{0.02in} {\bf Input:}
	loss function $f(\mathbf{w})$ with parameters $\mathbf{w}$; learning rate $\{\eta_t\}^T_{t=1}$; $\beta_2 = 0.999$\\
	\hspace*{0.02in} {\bf Output:} trained parameter
	\begin{algorithmic}[1]
		\State $\mathbf{m}_0 \gets \mathbf{0}$\enspace /* Initialize first-order momentum */
		\State $\mathbf{v}_0 \gets \mathbf{0}$, $\mathbf{\hat{v}}_0 \gets \mathbf{0}$\enspace /* Initialize second-order momentum */
		\For {$t = 1,2,\dots , T$}
		\State $\mathbf{g}_t = \nabla f_t(\mathbf{w}_t)$,
		\vspace{0.02in}
		\State $\theta = \left\langle\frac{\mathbf{m}_{t-1}}{\sqrt{\mathbf{\hat{v}}_{t-1}}}, \mathbf{g}_t\right\rangle$ \enspace /* The operator $\left\langle\cdot ,\cdot \right\rangle$ represents the angle between two vectors. */
		\If {$\theta \in [0,\frac{\pi}{2})$}
		\State $\beta_{1,t} = \sin\theta / K + \epsilon$
		\Else
		\State $\beta_{1,t} = 1 / K$ \enspace /* Here, $K = \frac{10(2+\pi)}{2\pi}$. */
		\EndIf

		\vspace{0.02in}
		\State $\mathbf{m}_t = (1-\beta_{1,t})\cdot\mathbf{m}_{t-1} + \beta_{1,t}\cdot\mathbf{g}_t$
		\vspace{0.02in}
		\State $\mathbf{v}_t = \beta_2\cdot\mathbf{v}_{t-1} + (1-\beta_2)\cdot\mathbf{g}_t\odot\mathbf{g}_t$\enspace /* $\odot$ is element-wise multiplication.*/
		\vspace{0.02in}
		\State $\mathbf{\hat{v}}_t = \max\{\mathbf{\hat{v}}_{t-1}, \mathbf{v}_t\}$
		\State $d_t = \min\{0.5\cos\theta, 0\}$
		\vspace{0.02in}
		\State $\mathbf{\Delta}_t = d_t\cdot\mathbf{\Delta}_{t-1} - \eta_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}}$
		\State $\mathbf{w}_t = \mathbf{w}_{t-1} + \mathbf{\Delta}_t$
		\EndFor

		\State \Return $\mathbf{w}_T$
	\end{algorithmic}\label{alg:cwm}

\end{algorithm}
{\text{DEAM}} introduces two new terms in the learning process: (1) the adaptive momentum weight $\beta_{1,t}
, and (2) the``backtrack term" $d_{t}$. In the $t_{th}$ training iteration, both $\beta_{1,t}$ and $d_t$ are calculated based on the ``discriminative angle" $\theta$, which is the angle between previous $\mathbf{m}_{t-1}/\sqrt{\mathbf{\hat{v}}_{t-1}}$ and current gradient $\mathbf{g}_t$ (since essentially both $\mathbf{m}_{t-1}/\sqrt{\mathbf{\hat{v}}_{t-1}}$ and $\mathbf{g}_t$ are vectors, there exists an angle between them). Here, $\mathbf{m}$ is the first-order momentum that records the exponential moving average of historical gradients; $\mathbf{v}$ is the exponential moving average of the squared gradients, which is called the second-order momentum. In the following parts of this paper, we will denote $\mathbf{m}_{t-1}/\sqrt{\mathbf{\hat{v}}_{t-1}}$ as the ``update volume'' in the ${(t-1)}_{th}$ iteration.
Formally, $\beta_{1,t}$ determines the weights of previous first-order momentum $\mathbf{m}_{t-1}$ and current gradient $\mathbf{g}_t$ when calculation the present $\mathbf{m}_t$. Meanwhile, the backtrack term $d_t$ represents the returning step of the previous update on parameters. We can notice that in each iteration, after the $\theta$ has been calculated, the $\beta_{1,t}$ and $d_t$ are directly obtained according to the $\theta$.
In this way, we can calculate appropriate $\beta_{1,t}$ as the discriminative angle changes. The $d_t$ term balances between the historical update term $\mathbf{\Delta}_{t-1}$ (defiend in Algorithm~\ref{alg:cwm}) and the current update volume $\mathbf{m}_t /\sqrt{\mathbf{\hat{v}}_t}$ when computing $\mathbf{\Delta}_t$. In the proposed {\text{DEAM}}, $\beta_{1,t}$ and $d_t$ terms can collaborate with each other and achieve faster convergence.

\subsection{Adaptive Momentum Weight $\beta_{1,t}$
\subsubsection{Motivation}\label{subsubsec: beta_motivation}
In the {\textsc{ADAM}}~\cite{ADAM} paper, (the first-order) momentum's weight (i.e., $\beta_1$) is a pre-specified fixed value, and commonly $\beta_1 = 0.9$. It has been used in many applications and the performance can usually meet the expectations. However, this setting is not applicable in some situations. For example, for the case
\begin{equation}
f(x,y) = x^2 + 4y^2 ,
\end{equation}
where $x$ and $y$ are two variables, it is obvious that $f$ is a convex function. If $f(x,y)$ is the objective function to optimize, we try to use {\textsc{ADAM}} to find its global optima.
\begin{figure}[t]
	\vspace{-20pt}
	\centering
	\includegraphics[width=6cm,height=3.75cm]{./counter-example2.pdf}
	\caption{The update routes of {\textsc{ADAM}} with $\beta_1 = 0.9$ (the blue line) and $\beta_1 = 0.0$ (the red line).}
	\label{fig:counter-example}
\end{figure}

\begin{figure*}[t]
	\centering
	\subfigure[The "zig-zag" route]{\label{fig:zigzag}
		\begin{minipage}[l]{0.6\columnwidth}
			\centering
			\includegraphics[width=0.8\textwidth,height=3cm]{./backtrack_11.pdf}
		\end{minipage}
	}

	\subfigure[Axis decomposition]{\label{fig:axis_decomp}
		\begin{minipage}[l]{0.6\columnwidth}
			\centering
			\includegraphics[width=0.8\textwidth,height=3cm]{./backtrack_2.pdf}
		\end{minipage}
	}

	\subfigure[Example: When $\theta = 180^\circ$]{\label{fig:theta180}
		\begin{minipage}[l]{0.6\columnwidth}
			\centering
			\includegraphics[width=0.8\textwidth,height=3cm]{./theta180_2.pdf}
		\end{minipage}
	}
	\caption{Some examples about $d_t$}
\end{figure*}
Let's assume ADAM starts the variable search from $(-4,-1)$ (i.e., the initial variable vector is $\mathbf{w}_0 = (-4,-1)^\top$) and the initial learning rate is $\eta_1 = 1$. Different choices of $\beta_1$ will lead to very different performance of ADAM. For instance, in Figure~\ref{fig:counter-example}, we illustrate the update routes of ADAM with $\beta_1 = 0.9$ and $\beta_1 = 0.0$ as the blue and red lines, respectively. In Figure~\ref{fig:counter-example}, the ellipse lines are the contour lines of $f(x,y)$, and points on the same line share the same function value. We can observe that after the first updating, both of the two approaches will update variables to $(-3, 0)$ point (i.e., the updated variable vector will be $\mathbf{w}_1 = (-3, 0)^\top$). In the second step, since the current gradient $\mathbf{g}_2 = (-6, 0)^\top$, the {\textsc{ADAM}} with $\beta_1 = 0.0$ will update variables in the $(1, 0)$ direction. Meanwhile, for the {\textsc{ADAM}} with $\beta_1 = 0.9$, its $\mathbf{m}_2 $ is computed by integrating $\mathbf{m}_1$ and $\mathbf{g}_2$ together (whose weights are $\beta_1$ and $1-\beta_1$, respectively). Therefore the updating direction of it will be more inclined to the previous direction instead. Compared with ADAM with $\beta_1 = 0.0$, the ADAM with $\beta_1 = 0.9$ takes much more iterations until converging.

From the analysis above, we can observe that, a careful tuning and updating of $\beta_1$ in the learning process can be crucial for the performance of {\textsc{ADAM}}. However, by this context so far, there still exist no effective approaches for guiding the parameter tuning yet. To deal with this problem, {\text{DEAM}} introduces the concept of discriminative angle $\theta$ for computing $\beta_1$ automatically as follows.

\subsubsection{Mechanism}
The momentum weight $\beta_1$ will be updated in each iteration in {\text{DEAM}}, and we can denote its value computed in the $t_{th}$ iteration as $\beta_{1,t}$ formally.
Essentially, in the $t_{th}$ iteration of the training process, both the previous update volume
and $\mathbf{g}_t$ are vectors (or directions), and these directions directly decide the updating process. Thus we try to extract their relation with the help of angle, and subsequently determine the weight $\beta_{1,t}$ (or $1 - \beta_{1,t}$) by the angle.

In Algorithm~\ref{alg:cwm}, the discriminative angle $\theta$ in the $t_{th}$ iteration is calculated by
\begin{small}
	\begin{equation}
	\theta = \left\langle -\frac{\mathbf{m}_{t-1}}{\sqrt{\mathbf{\hat{v}}_{t-1}}}, -\mathbf{g}_t \right\rangle = \left\langle\frac{\mathbf{m}_{t-1}}{\sqrt{\mathbf{\hat{v}}_{t-1}}}, \mathbf{g}_t\right\rangle
	\end{equation}
\end{small}
Here, the operator $\left\langle\cdot, \cdot\right\rangle$ denotes the angle between two vectors. This expression is easy to understand, since the $-\mathbf{m}_{t-1} / \sqrt{\mathbf{\hat{v}}_{t-1}}$ can represent the updating direction of $(t-1)_{th}$ iteration in {\text{AMSGrad}}, meanwhile $-\mathbf{g}_t$ is the reverse of the present gradient. So we can simplify it as $\theta = <\frac{\mathbf{m}_{t-1}}{\sqrt{\mathbf{\hat{v}}_{t-1}}}, \mathbf{g}_t>$. If $\theta$ is close to zero (denoted by $\theta\rightarrow 0^\circ$), the $\mathbf{m}_{t-1} / \sqrt{\mathbf{\hat{v}}_{t-1}}$ (previous update volume) and $\mathbf{g}_t$ are almost in the same direction, and the weights for them will not be very important. Meanwhile, if $\theta$ approaches $180^\circ$ (denoted by $\theta\rightarrow 180^\circ$), the previous update volume and $\mathbf{g}_t$ will be in totally reverse directions. This means in the current step, the previous momentum term is already in a wrong direction. Therefore, to rectify this error of the last momentum, {\text{DEAM}} proposes to assign the current gradient's weight (i.e., $\beta_{1,t}$ in our paper) with a larger value instead. As the $\beta_{1,t}$ varies when $\theta$ changes from $0^\circ$ to $180^\circ$, we intend to define $\beta_{1,t}$ with the following equation:
\begin{small}
	\begin{equation}
	\beta_{1,t} =  \begin{cases}
	\sin \ K + \epsilon &  \mbox{$\theta \in [0, \frac{\pi}{2})$} \\
	1 / K & \mbox{$\theta\in [\frac{\pi}{2}, \pi]$}
	\end{cases}
	\end{equation}\label{equ:beta_1}
\end{small}
where $K = 10(2+\pi)/2\pi$ and $\epsilon$ is a very small value (e.g., $\epsilon = 0.001$).
In the equation above, the threshold of the piecewise function is $\theta = \pi/2$, because $\sin\theta$ comes to the maximum at this point and goes down when $\theta > \frac{\pi}{2}$. If $\frac{\pi}{2}\leq \theta\leq \pi$, which is exactly the situation $\theta\rightarrow 180^\circ$ we discussed above, we intend to keep $\beta_{1,t}$ in a relatively large value. The reason we rescale $\sin\theta$ by ${1}/{K}$ is that directly applying $\beta_{1,t}= \sin\theta$ will overweight $\mathbf{g}_t$, which may cause fluctuations on the update routes. The value of $K$ is determined by:
\begin{small}
	\begin{equation}
	K = 10(\int_{0}^{\frac{\pi}{2}} \sin\theta d\theta + \int_{\frac{\pi}{2}}^{\pi} 1 d\theta) = \frac{10(2+\pi)}{2\pi}
	\end{equation}\label{equ:K}
\end{small}
In the equation above, assume $\theta$ is randomly distributed on $[0, \pi]$, in this calculation we can get
\begin{small}
	\begin{equation}
	\mathbb{E}[\beta_{1,t}] = \frac{1}{\pi}\int_{0}^{\pi} \beta_{1,t}(\theta) d\theta = 0.1
	\end{equation}\label{equ:E[beta_1]}
\end{small}
In other words, the expectation of $\beta_{1,t}$ (i.e., $\mathbb{E}(\beta_{1,t})$) will be identical to the $\beta_1$ used in {\textsc{ADAM}}~\cite{ADAM}.
After obtaining $\beta_{1,t}$, it will be applied to calculating $\mathbf{m}_t$ as shown in Algorithm~\ref{alg:cwm}. In this way, we have achieved momentum with adaptive weights, and this weight is automatically computed during the training process, fewer hyperparameters will be involved.

\subsection{Backtrack Term $d_t$}
To further speed up the convergence rate, we employ a novel backtrack mechanism for {\text{DEAM}}. As a mechanism computed based on the discriminative angle $\theta$, the backtrack term allows {\text{DEAM}} to eliminate redundant update in each iteration.
Besides, according to our following analysis, the backtrack term $d_t$ virtually collaborates with the $\beta_{1,t}$ term to further accelerate the convergence of the training process.
\vspace{-5pt}
\subsubsection{Motivation}\label{sec:motivation}
When optimizer (e.g., {\textsc{ADAM}}) updates variables of the loss function (e.g., $f(x, y)$), some update routes will look like the black arrow lines shown in Figure~\ref{fig:zigzag}, especially when the discriminative angle $\theta$ is larger than $90^\circ$. We call this phenomenon the "zig-zag" route. In Figure~\ref{fig:zigzag}, it shows the update routes of a 2-dimension function. Each black arrow line in the figure represents the variables' update in each epoch; the red dashed line is the direction of the update routes; the $\theta$ is the discriminative angle. If $\theta \geq 90^\circ$, the "zig-zag" phenomenon will appear severely, which may lead to slower convergence speed. The main reason is when $\theta \geq 90^\circ$, if we map two neighboring update directions onto the coordinate axes, there will be at least one axis of the directions being opposite. This situation is shown in Figure~\ref{fig:axis_decomp}. For the example of a function with 2-dimension variables, the update volume $\mathbf{m}_1/\sqrt{\mathbf{v}_1}$ can be decomposed into $(x_1, y_1)^\top$ in Figure~\ref{fig:axis_decomp}, and the same with $\mathbf{m}_2/\sqrt{\mathbf{v}_2}$. We can notice that $y_1$ and $y_2$ are in the opposite directions, so the first and second steps practically have inverse updates subject to the $y$ axis. We attribute this situation to the over update (or redundant update) of the first step. Therefore the backtrack term $d_t$ is proposed to restrict this situation.

\subsubsection{Mechanism}\label{subsubsec:dt_m}
Since the redundant update situation is caused by over updating of the previous iteration, simply we intend to deal with it through a backward step. Meanwhile, during the updating process of variables, not every step will suffer from the redundant update: if $\theta\rightarrow 0^\circ$, the updating process becomes smooth, not like the situation shown in Figure~\ref{fig:zigzag}. Besides, from the analysis above we conclude that if $\theta \geq 90^\circ$, there will be at least one dimension involves the redundant update. Thus, in the $t_{th}$ iteration we quantify $d_t$ as the following equation:
\begin{small}
	\begin{equation}
	d_t = \min\{0.5\cos\theta, 0\}
	\end{equation}\label{equ:d_t}
\end{small}
and we rewrite the updating term with backtrack in {\text{DEAM}} as
\begin{small}
	\begin{equation}
	\mathbf{\Delta}_t = d_t\cdot\mathbf{\Delta}_{t-1} - \eta_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}}
	\end{equation}\label{equ:Delta_t}
\end{small}
where $\theta$ is the discriminative angle and $\mathbf{\Delta}_t$ is the updating term in Algorithm~\ref{alg:cwm}. By designing $d_t$ in this way, when $\theta\rightarrow 0^\circ$, $d_t = 0$ and there is no backward step, the updating term $\mathbf{\Delta}_t = -\eta_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}}$ is similar to {\text{AMSGrad}}; when $\theta\rightarrow 180^\circ$, $d_t = 0.5\cos\theta$ and comes to the maximum value when $\theta = 180^\circ$. The reason that $\cos\theta$ is rescaled by $0.5$ is that: in Figure~\ref{fig:theta180}, $\mathbf{w}_{t-1}$ and $\mathbf{w}_t$ are the variables updated by {\text{DEAM}} without $d_t$ term in the $(t-1)_{th}$ and $t_{th}$ iterations respectively. If the backtrack mechanism is implemented, in the $(t+1)_{th}$ iteration, since $\theta = 180^\circ$, firstly $d_t = 0.5\cos\theta \to -0.5$ makes the backtrack to the $\mathbf{w'}_t$ point (the middle point of $\mathbf{w}_{t-1}$ and $\mathbf{w}_t$). Thus, this backtrack step allows the variable to further approach the optima.

By implementing the backtrack term $d_t$, {\text{DEAM}} can combine it with the adaptive momentum weight $\beta_{1,t}$ to achieve the collaborating of them. For the situation of large discriminative angle ($\theta > 90^\circ$), both $\beta_{1,t}$ and $d_t$ in the current step can make corrections to the last update. Since when $\theta > 90^\circ$, the last update is in conflict direction compared with the current gradient, and  $\beta_{1,t}$ will increase to allocate a large weight for the present gradient, which subsequently corrects the previous step. Meanwhile, the $d_t$ will also conduct a backward step of to further rectify the last update.

\begin{figure*}[t]
	\centering
	\subfigure[Train loss on ORL]{
		\begin{minipage}[l]{0.48\columnwidth}
			\centering
			\includegraphics[width=0.85\textwidth]{./lr_orl_train_loss_big.pdf}

		\end{minipage}
		\label{fig:orl_1}
	}
	\subfigure[Test loss on ORL]{
		\begin{minipage}[l]{0.48\columnwidth}
			\centering
			\includegraphics[width=0.85\textwidth]{./lr_orl_test_loss_big.pdf}

		\end{minipage}
		\label{fig:orl_2}
	}
	\subfigure[Train loss on MNIST]{ \label{fig:mnist_train}
		\begin{minipage}[l]{0.48\columnwidth}
			\centering
			\includegraphics[width=0.85\textwidth]{./dnn_mnist_train_loss_big.pdf}

		\end{minipage}
		\label{fig:dnn_1}
	}

	\subfigure[Test loss on MNIST]{ \label{fig:mnist_test}
		\begin{minipage}[l]{0.48\columnwidth}
			\centering
			\includegraphics[width=0.85\textwidth]{./dnn_mnist_test_loss_big.pdf}

		\end{minipage}
		\label{fig:dnn_2}
	}
	\caption{Results of Logistic Regression and DNN structures}\label{fig:dnn}
\end{figure*}
\subsection{Theoretical Analysis}\label{sec: convergence_analysis}
In this part, we give the detailed analysis on the convergence of our  {\text{DEAM}} algorithm. According to~\cite{ADAM,amsgrad,online_convex_program,adafom}, given an arbitrary sequence of convex objective functions $f_1(\mathbf{w}), f_2(\mathbf{w}),\dots, f_T(\mathbf{w})$, we intend to evaluate our algorithm using the regret function, which is denoted as:
\begin{small}
	\begin{equation}
	R(T) = \sum_{t=1}^{T}[f_t(\mathbf{w}_t) - f_t(\mathbf{w}^*)]
	\end{equation}\label{equ:Delta_t}
\end{small}
where $\mathbf{w}^*$ is the globally optimal point. In the following Theorem~\ref{thm:converge}, we will show that the above regret function is bounded. Before proving the Theorem~\ref{thm:converge}, there are some definitions and lemmas as the pre-requisites.
\begin{defn} \label{def:convex}

	If a function $f: \mathbb{R}^d \rightarrow \mathbb{R}$ is convex, then $\forall x, y\in\mathbb{R}^d$, $\forall \phi\in[0,1]$, we have
	\[
	f(\phi x + (1-\phi)y)\leq\phi f(x) + (1-\phi)f(y)
	\]
\end{defn}
\begin{defn}\label{lemma:bound}
	If a function $f: \mathbb{R}^d\rightarrow \mathbb{R}$ is convex, then $\forall x, y\in\mathbb{R}^d$ we have
	\[
	f(y)\geq f(x) + \nabla f(x)^\top (y-x)
	\]
\end{defn}

\begin{lemma}\label{lemma:m_t}
	Assume that the function $f_t$ has bounded gradients, $\left\|\nabla f_t(\mathbf{w})\right\|_\infty\leq G_\infty$. Let  $\mathbf{m}_{t, i}$ represents the $i_{th}$ element of $\mathbf{m}_t$ in {\text{DEAM}}, then the $\mathbf{m}_{t, i}$ is bounded by
	\[
	\mathbf{m}_{t,i}\leq \frac{(1-\epsilon_0)G_\infty}{K(1-\lambda)}
	\]
\end{lemma}
\begin{proof}
	Let $g_t = \nabla f_t(\mathbf{w})$.
	According to the definition of $\mathbf{m}_{t,i}$ in our algorithm,
	{\scriptsize
		\begin{align*}
		\mathbf{\mathbf{m}}_{t,i} &= \sum_{j=1}^{t}\beta_{1,j}\prod_{l=1}^{t-j}(1-\beta_{1,t-l+1})g_{j,i}\\

		&\leq \frac{G_\infty}{K}\sum_{j=1}^{t}\prod_{l=1}^{t-j}(1-\epsilon)
		\leq \frac{G_\infty}{K}\sum_{j=1}^{t}(1-\epsilon_0)\lambda^{t-j}\leq \frac{(1-\epsilon_0)G_\infty}{K(1-\lambda)}
		\end{align*}
	}
	where $K$ and $\epsilon$ are the terms in Algorithm~\ref{alg:cwm}.
\end{proof}
For the following proof, $\mathbf{g}_t := \nabla f_t(\mathbf{w}_t)$ and $\mathbf{g}_{t, i}$ will represent the $i_{th}$ element of $\mathbf{g}_t\in\mathbb{R}^d$, and $\mathbf{g}_{1:t, i} = [\mathbf{g}_{1,i},\mathbf{g}_{2,i},\dots ,\mathbf{g}_{t,i}]$.
\begin{thm}\label{thm:converge}
	Assume $\{f_t\}_{t=1}^T$ have bounded gradients $\left\|\nabla f_t(\mathbf{w})\right\|_\infty\leq G_\infty$ for all $\mathbf{w}\in \mathbb{R}^d$, all variables are bounded by $\left\|\mathbf{w}_p - \mathbf{w}_q\right\|_2\leq D$ and $\left\|\mathbf{w}_p - \mathbf{w}_q\right\|_\infty\leq D_\infty$, $\forall p,q \in \{1,2,\dots,T\}$, $\eta_t = \eta/\sqrt{t}$, $\gamma_1 = (1-\epsilon_0) / \sqrt{\beta_2}$ and satisfies $\gamma_1<1$, $\epsilon = 1-(1-\epsilon_0)\lambda^{t-1}, \lambda\in(0,1)$. Our proposed algorithm can achieve the following bound on regret:
	{\scriptsize
		\begin{align*}

		R(t)&\leq \frac{D^2}{\epsilon_0\eta}\sum_{i=1}^{d}\sqrt{T\mathbf{\hat{v}}_{T,i}}
		+ \frac{(1-\epsilon_0)^2G_\infty D_\infty d}{K(1-\lambda)^2\epsilon_0}\\
		&+ \frac{\eta \sqrt{1+\log T}}{2\epsilon^2_0(1-\gamma_1)\sqrt{1-\beta_2}}\sum_{i=1}^{d}\left\|\mathbf{g}_{1:T,i}\right\|_2
		\end{align*}
	}

\end{thm}
\begin{proof}
	According to Definition~\ref{lemma:bound}, for $\forall t \in\{1,2,\dots,T\}$, we have
	{\scriptsize
		\begin{align*}
		f_t(\mathbf{w}_t) - f_t(\mathbf{w}^*) \leq\nabla f_t(\mathbf{w}_t)^\top (\mathbf{w}_t - \mathbf{w}^*)

		= \sum_{i=1}^{d} \mathbf{g}_{t, i}(\mathbf{w}_{t,i} - \mathbf{w}_{i}^*)
		\end{align*}
	}

	From the definition of $\mathbf{\Delta}_{t}$ in the updating rule of {\text{DEAM}}, we know it is equal to multiplying the learning rate $\eta_t$ in some iterations by a number in $[0.5, 1]$, which means $\mathbf{w}_{t+1} = \mathbf{w}_t - \hat{\eta}_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}};\hat{\eta}_t = \mu_t\cdot\eta_t$, where $\mu_t\in [0.5, 1]$.
	If we first focus on the $i_{th}$ element of $\mathbf{w}_t$, we can get
	{\scriptsize
		\begin{align*}
		&(\mathbf{w}_{t+1,i} - \mathbf{w}^*_i)^2= (\mathbf{w}_{t,i} - \mathbf{w}^*_i - \hat{\eta}_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}})^2 =(\mathbf{w}_{t,i} - \mathbf{w}^*_i)^2 \\

		-& 2\hat{\eta}_t(\frac{(1 - \beta_{1,t})}{\sqrt{\mathbf{\hat{v}}_{t,i}}}\mathbf{m}_{t-1,i} + \frac{\beta_{1,t}}{\sqrt{\mathbf{\hat{v}}_{t,i}}}\mathbf{g}_{t,i})\cdot(\mathbf{w}_{t,i} - \mathbf{w}^*_i)+ (\hat{\eta}_t\cdot\frac{\mathbf{m}_t}{\sqrt{\mathbf{\hat{v}}_t}})^2
		\end{align*}
	}
	Then,
	{\scriptsize
		\begin{align*}
		&2\hat{\eta}_t\cdot\frac{\beta_{1,t}}{\sqrt{\mathbf{\hat{v}}_{t,i}}}\mathbf{g}_{t,i} (\mathbf{w}_{t,i} - \mathbf{w}^*_i) = (\mathbf{w}_{t,i} -\mathbf{w}^*_i)^2 - (\mathbf{w}_{t+1,i} -\mathbf{w}_i^*)^2\\
		&- 2\hat{\eta}_t\cdot\frac{(1 - \beta_{1,t})}{\sqrt{\mathbf{\hat{v}}_{t,i}}}\mathbf{m}_{t-1,i}\cdot(\mathbf{w}_{t,i} - \mathbf{w}^*_i) + \hat{\eta}^2_t\cdot\frac{\mathbf{m}^2_{t,i}}{\mathbf{\hat{v}}_{t,i}}
		\end{align*}
	}
	So we can obtain
	{\scriptsize
		\begin{align}
		\mathbf{g}_{t,i} (\mathbf{w}_{t,i} - \mathbf{w}^*_i) &= \frac{\sqrt{\mathbf{\hat{v}}_{t,i}}}{2\hat{\eta}_t\beta_{1,t}} [(\mathbf{w}_{t,i} - \mathbf{w}^*_i)^2 - (\mathbf{w}_{t+1,i} - \mathbf{w}^*_i)^2]\\
		&~~~~ + \frac{(1-\beta_{1,t})}{\beta_{1,t}}\mathbf{m}_{t-1,i}(\mathbf{w}^*_i - \mathbf{w}_{t,i})\\
		&~~~~ +\frac{\hat{\eta}_t}{2\beta_{1,t}} \cdot\frac{\mathbf{m}^2_{t,i}}{\sqrt{\mathbf{\hat{v}}_{t,i}}}
		\end{align}\label{fmla:1}
	}
	For the right part of $(9)$ in the above formula, if we sum it from $t=1$ to $t = T$,
	{\scriptsize
		\begin{align*}
		&\sum_{t=1}^{T}\frac{\sqrt{\mathbf{\hat{v}}_{t,i}}}{2\hat{\eta}_t\beta_{1,t}}[(\mathbf{w}_{t,i} - \mathbf{w}^*_i)^2 - (\mathbf{w}_{t+1,i} - \mathbf{w}^*_i)^2]\\

		\leq& \frac{1}{\epsilon_0}\{(\mathbf{w}_{1,i} - \mathbf{w}^*_i)^2\cdot\frac{\sqrt{\mathbf{\hat{v}}_{1,i}}}{\eta_1}

		+\dots + (\mathbf{w}_{T,i} - \mathbf{w}^*_i)^2(\frac{\sqrt{\mathbf{\hat{v}}_{T,i}}}{\eta_T} - \frac{\sqrt{\mathbf{\hat{v}}_{T-1,i}}}{\eta_{T-1}})\}\\
		\leq& \frac{D^2}{\epsilon_0\eta}\sqrt{T\mathbf{\hat{v}}_{T,i}}
		\end{align*}
	}
	The first inequality is satisfied because of the line 13 in Algorithm~\ref{alg:cwm}. For the $(10)$ in the formula, if we sum it from $t=1$ to $t = T$,
	{\scriptsize
		\begin{align*}
		&\sum_{t=1}^{T}\frac{(1-\beta_{1,t})}{\beta_{1,t}}\mathbf{m}_{t-1,i}(\mathbf{w}^*_i - \mathbf{w}_{t,i})
		\leq\frac{(1-\epsilon_0)G_\infty D_\infty}{K(1-\lambda)\epsilon_0} \sum_{t=1}^{T}(1-\beta_{1,t})\\
		\leq& \frac{(1-\epsilon_0)G_\infty D_\infty}{K(1-\lambda)\epsilon_0}\sum_{t=1}^{T}(1 - \epsilon)
		= \frac{(1-\epsilon_0)G_\infty D_\infty}{K(1-\lambda)\epsilon_0}\sum_{t=1}^{T}(1 - \epsilon_0)\lambda^{t-1}\\
		\leq& \frac{(1-\epsilon_0)^2G_\infty D_\infty}{K(1-\lambda)^2\epsilon_0}
		\end{align*}
	}

	\begin{figure*}[t]

		\centering
		\subfigure[Train loss on ORL]{\label{fig:orl_train}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_orl_train_loss_big.pdf}

			\end{minipage}
		}

		\subfigure[Test loss on ORL]{ \label{fig:orl_test}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_orl_test_loss_big.pdf}

			\end{minipage}
		}
		\vspace{-5pt}
		\subfigure[Train loss on MNIST]{ \label{fig:mnist_cnn_train}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_mnist_train_loss_big.pdf}

			\end{minipage}
		}

		\subfigure[Test loss on MNIST]{ \label{fig:mnist_cnn_test}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_mnist_test_loss_big.pdf}

			\end{minipage}
		}

		\subfigure[Train loss on CIFAR]{ \label{fig:cifar_train}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_cifar_train_loss_big.pdf}

			\end{minipage}
		}

		\subfigure[Test loss on CIFAR]{ \label{fig:cifar_test}
			\begin{minipage}[l]{0.65\columnwidth}
				\centering
				\includegraphics[width=0.85\textwidth]{./cnn_cifar_test_loss_big.pdf}

			\end{minipage}
		}
		\vspace{-10pt}
		\caption{Results of CNN structure}\label{fig:cnn}

	\end{figure*}

	\begin{table*}[t]
		\renewcommand\arraystretch{1.1}
		\scriptsize
		\centering
		\begin{threeparttable}
			\caption{Running time of {\text{DEAM}} and comparison methods (the unit of values is second)}
			\begin{tabular}{|c|c|c|c|c|c|}
				\hline
				\multirow{2}*{\tabincell{c}{\textbf{Comparison}\\\textbf{Methods}}}&\multicolumn{5}{|c|}{Running time on all models}\\

				\cline{2-6}
				&\tabincell{c}{Logistic Regression \\ on ORL}&DNN on MNIST&CNN on ORL&CNN on MNIST&CNN on CIFAR-10\\
				\cline{1-6}
				{\text{DEAM}}&\textbf{38}&\textbf{302}&\textbf{35064}&\textbf{11679}&\textbf{57761}\\
				\hline
				{\textsc{ADAM}}&102&664&47418&21775&67584\\
				\hline
				{\text{RMSProp}}&48&307&36722&11997&84305\\
				\hline
				{\text{AdaGrad}}&$>200$&667&$>100000$&$>50000$&$>100000$\\
				\hline
				SGD&$>200$&346&$>100000$&16985&67564\\
				\hline
			\end{tabular}

			\label{tab:running_time}
		\end{threeparttable}
	\end{table*}

	The first inequality is according to Lemma~\ref{lemma:m_t}.	Finally, we will infer the $(11)$ in previous formula. According to the Lemma 2 of~\cite{amsgrad}, we have
	{\scriptsize
		\begin{align*}
		\small
		&\sum_{t=1}^{T}\frac{\hat{\eta}_t}{2\beta_{1,t}}\cdot\frac{\mathbf{m}^2_{t,i}}{\sqrt{\mathbf{\hat{v}}_{t,i}}}\leq \frac{1}{2\epsilon_0}\sum_{t=1}^{T}\eta_t\frac{\mathbf{m}^2_{t,i}}{\sqrt{\mathbf{v}_{t,i}}}
		\leq \frac{\eta}{2\epsilon_0}\sum_{t=1}^{T}\frac{1}{\sqrt{t}}\cdot\frac{\mathbf{m}^2_{t,i}}{\sqrt{\mathbf{v}_{t,i}}}\\

		\leq& \frac{\eta}{2\epsilon_0} \sum_{t=1}^{T}\frac{(\sum_{j=1}^{t}(1-\epsilon_0)^{t-j})(\sum_{j=1}^{t}(1-\epsilon_0)^{t-j}\mathbf{g}^2_{j,i})}{\sqrt{t((1-\beta_2)\sum_{j=1}^{t}\beta_2^{t-j}\mathbf{g}^2_{j,i})}}\\

		\leq& \frac{\eta }{2\epsilon^2_0\sqrt{1-\beta_2}}\sum_{t=1}^{T}\frac{1}{\sqrt{t}}\sum_{j=1}^{t}\frac{(1-\epsilon_0)^{t-j}\mathbf{g}^2_{j,i}}{\sqrt{\beta_2^{t-j}\mathbf{g}^2_{j,i}}}\\

		\leq& \frac{\eta }{2\epsilon^2_0\sqrt{1-\beta_2}}\sum_{t=1}^{T}|\mathbf{g}_{t,i}|\sum_{j=t}^{T}\frac{\gamma_1^{j-t}}{\sqrt{t}}

		\leq \frac{\eta \sqrt{1+\log T}}{2\epsilon^2_0(1-\gamma_1)\sqrt{1-\beta_2}}\left\|\mathbf{g}_{1:T,i}\right\|_2
		\end{align*}
	}
	In the above inequalities, some inferences are based on Cauchy-Schwarz Inequality. Therefore, the final bound of $R(T)$ can be expressed as
	{\scriptsize
		\begin{align*}
		R(T) &\leq \frac{D^2}{\epsilon_0\eta}\sum_{i=1}^{d}\sqrt{T\mathbf{\hat{v}}_{T,i}} + \frac{(1-\epsilon_0)^2G_\infty D_\infty d}{K(1-\lambda)^2\epsilon_0}\\ &+ \frac{\eta \sqrt{1+\log T}}{2\epsilon^2_0(1-\gamma_1)\sqrt{1-\beta_2}}\sum_{i=1}^{d}\left\|\mathbf{g}_{1:T,i}\right\|_2
		\end{align*}
	}
\end{proof}
For the bound term, as $T\rightarrow +\infty$, $\frac{R(T)}{T}\rightarrow 0$ and we can infer that $\displaystyle\lim_{T\rightarrow \infty} [f_t(\mathbf{w}_t) - f_t(\mathbf{w}^*)] = 0$, which means the proposed algorithm can finally converge.

\section{Experiments}\label{sec:experiment}
We have applied the {\text{DEAM}} algorithm on multiple popular machine learning and deep learning structures, both convex and non-convex situations. To show the advantages of the algorithm, we compare it with various popular optimization algorithms, including {\textsc{ADAM}}~\cite{ADAM}, RMSProp~\cite{rmsprop}, AdaGrad~\cite{adagrad} and SGD. For all the experiments, the loss function (objective function) we have selected is the cross-entropy loss function, and the size of the minibatch is 128. Besides, the learning rate is 0.0001.

\subsection{Experiment Settings and Results}

\noindent\textbf{Logistic Regression}: We firstly evaluate our algorithm on the multi-class logistic regression model, since it is widely used and owns a convex objective function. We conduct logistic regression on the ORL dataset~\cite{ORL}. ORL dataset consists of face images of 40 people, each person has ten images and each image is in the size of $112\times 92$. The loss of objective functions on both training set and testing set are shown in Figure~\ref{fig:orl_1}, \ref{fig:orl_2}.

\noindent\textbf{Deep Neural Network}: We use deep neural network (DNN) with two fully connected layers of 1,000 hidden units and the Relu~\cite{relu} activation function. The dataset we use is MNIST~\cite{GBDR}. The MNIST dataset includes 60,000 training samples and 10,000 testing samples, where each sample is a $28\times 28$ image of hand-written numbers from 0 to 9. Result are exhibited in Figure~\ref{fig:dnn_1}, \ref{fig:dnn_2}.

\noindent\textbf{Convolutional Neural Network}: The CNN models in our experiments are based on the LeNet-5~\cite{GBDR}, and it is implemented on multiple datasets: ORL, MNIST and CIFAR-10~\cite{CIFAR}. The CIFAR-10 dataset consists of 60,000 $32\times 32$ images in 10 classes, with 6,000 images per class. For different datatsets, the structures of CNN models are modified: for the ORL dataset, the CNN model has two convolutional layers with 16 and 36 feature maps of $5\time 5$ kernels and 2 max-pooling layers, and a fully connected layer with 1024 neurons; for the MNIST dataset, the CNN structure follows the LeNet-5 structure in~\cite{GBDR}; for CIFAR-10 dataset, the CNN model consists of three convolutional layers with 64, 128, 256 kernels respectively, and a fully connected layer having 1024 neurons. All experiments apply Relu~\cite{relu} activation function. The results are shown in Figure~\ref{fig:cnn}.

We can observe that {\text{DEAM}} converges faster than other widely used optimization algorithms in all the cases. Within the same number of epoches, {\text{DEAM}} can converge to the lowest loss on both the training set and test set.

\subsection{Analysis of the Backtrack Mechanism}

To show the effectiveness of the backtrack term $d_t$, we also carry out the experiments of {\text{DEAM}} without $d_t$ term, and exhibit the results in Figures~\ref{fig:dnn}, \ref{fig:cnn}. The results indicate that after applying $d_t$ term, the converging speed will become slightly faster. To thoroughly prove the effectiveness of our proposed $d_t$ term, we also compare it with other definitions of backtrack terms e.g., $d_t = 0.5\cos\theta$, and exhibit the results in Figure~\ref{fig:d_t_other}. In Figure~\ref{fig:d_t_other}, $d_t = sigmoid \_ based$ represents that $d_t = -1/(1+ e^{-(\theta - \frac{1}{2}\pi)}) + \frac{1}{2}$, and $d_t = tanh \_ based$ means that $d_t = -2/(1 + e^{-2(x - \frac{1}{2}\pi)}) + 1$.
\begin{figure}[t]

	\centering
	\includegraphics[width=5cm,height=3.5cm]{./orl_cnn_dt.pdf}
	\caption{$d_t$ terms in other definitions}
	\vspace{-5pt}
	\label{fig:d_t_other}
\end{figure}
Due to the limited sapce, here we only exhibit the results on ORL dataset. From the results we can observe that the $d_t$ definnition in {\text{DEAM}} achieves the best congerging performance, which means our setting of $d_t$ in {\text{DEAM}} is effective.

\subsection{Time-consuming Analysis}

We have recorded the running time of {\text{DEAM}} and other comparison algorithms in every experiment, and list them in the Table~\ref{tab:running_time}. The running time shown in Table~\ref{tab:running_time} contains ``$>$'', which means the model still does not converge at the specific time. From the results we can observe that in all of our experiments, {\text{DEAM}} finally converges within the smallest mount of time. From the results in Figures~\ref{fig:dnn},\ref{fig:cnn} and Table~\ref{tab:running_time}, we can conclude that {\text{DEAM}} can converge not only in fewer epochs, but using less running time. The device we used is the Dell PowerEdge T630 Tower Server, with 80 cores 64-bit Intel Xeon CPU E5-2698 v4@2.2GHz. The total memory is 256 GB, with an extra (SSD) swap of 250 GB.

\section{Conclusion}\label{sec:conclusion}
In this paper, we have introduced a novel optimization algorithm, the {\text{DEAM}}, which implements the momentum with discriminative weights and the backtrack term. We have analyzed the advantages of the proposed algorithm and proved it by theoretical inference. Extensive experiments have shown that the proposed algorithm can converge faster than existing methods on both convex and non-convex situations, and the time consuming is better than existing methods. Not only the proposed algorithm can outperform other popular optimization algorithms, but fewer hyperparameters will be introduced, which makes the {\text{DEAM}} much more applicable.

\newpage
\bibliographystyle{named}


\section{Introduction}
	A qubit corresponds to the fundamental unit of quantum information in the same sense as the bit corresponds to the fundamental classical unit of information. When we are referring to quantum states with a Hilbert space with dimension $N > 2$, then we speak of \textit{qudits} \cite{rungta2001qudit, thew2004bell, de2002creating, thew2002qudit, neves2005generation, moreva2006realization, li2008generation}. The use of qudits opens new prospects for quantum information processing, such as the enhancement of security for quantum key distribution, the increase of channel capacity, higher noise resistance and quantum computing \cite{Ursin17, Torres12, kues2017chip, wang2018multidimensional, walborn2006quantum, wang2017qudit, wang2018proof, etcheverry2013quantum, gedik2015computational, niu2018qudit, chau2015quantum, islam2017provably, pal2018entanglement}.

	Manipulating states of qudits in interferometers at optical tables is a non-scalable and challenging task. Interferometers have been set for characterizing quantum states \cite{agnew2011tomography,Pimenta13,bent2015experimental,ikuta2017implementation,martinez2019experimental}, controlling, detecting or measuring entanglement of qudit systems \cite{kues2017chip,wang2018multidimensional, gutierrez2012experimental, gutierrez2014detection} or simulating noisy channels \cite{marques2015experimental,varga2018characterizing}. The use of photonic circuits written in glasses by femtosecond direct laser writing technique opens new possibilities for the generation, manipulation, and characterization of high dimensional states, besides the direct transmission of these states through an optical fiber \cite{Politi08,Marshall09,shadbolt2012generating,o2009photonic,caspani2017integrated,smith2009phase,Davis96,Streltsov01}. Using longitudinal interferometers, we can implement quantum operations in two qudits states, including transforming them from antisymmetric states to symmetric ones, making possible to construct interferometers for simulation of open systems and their characterization \cite{Walther12, Koch12,harris2017quantum,kerman2017multiloop}. Remote preparation of arbitrary qudits for their use in quantum communication protocols is another possible application for these circuits \cite{Fan08, Yin11,wang2015controlled,cao2016flexible}.

    Recently, the generation and analysis of a two qutrits system in a silicon photonic circuit generated with the lithographic technique was demonstrated \cite{Schaeff15}. In this domain, we can use such result to produce states of two qudits from an integrated photonic chip drawn on glass, after they receive photon pairs generated by spontaneous parametric down-conversion (SPDC). We aim here the construction of photonic circuits that, by implementing POVM  (Positive Operator-Valued Measure) elements, are able to perform minimum full quantum state tomography of qudits, defined as an experimental technique that allows the reconstruction of the density matrix of a quantum system in a unknown state with a minimum number of measurements \cite{Raymer94, Home2006, Riebe2006,six2016quantum,delaney2019quantum,xin2017quantum}. These techniques satisfy a demand in current computer and telecommunications technologies: the ability to execute quantum information protocols on small devices that integrate a larger hardware design. Our accomplishment is a contribution in the development of quantum information tasks in chips. It contributes to the technology needed to reach compact and portable devices for quantum protocols.

    Quantum state tomography on chip has been a subject of interest to various research groups. It was demonstrated, theoretically and experimentally, the accomplishment of quantum tomography in photonic chips through quantum walks \cite{Titchener16,titchener2018scalable}. The authors used in their experiment a circuit with N inputs and M outputs to characterize a photon-number state. The detection of the photons of all the outputs is done simultaneously, optimizing the experiment duration. Since it is not necessary to reconfigure the experimental apparatus, the occurrence of errors is also reduced. Circuits that do not require reconfiguration are called static circuits \cite{oren2017quantum}.

    In this work, we present an alternative form of performing full quantum tomography on photonic chips, by the implementation of POVM elements in static circuits, whose detection at all outputs is made simultaneously. All necessary POVM elements are implemented at the same time by the circuit and measured at the different photonic circuit exits simultaneously. All necessary measurement probabilities are obtained from the photon counts at the circuit outputs in one measurement time interval. This is important because it minimizes the noise introduced by the experimental apparatus during the detection time. The proposed tomography method uses SIC-POVM (Symmetric Informationally Complete Positive Operator-Valued Measure) elements, allowing a full quantum state tomography with the minimum number of measurements and with a symmetric POVM elements distribution in Hilbert space that optimizes the system density operator reconstruction \cite{Pimenta13,PaivaSanchez,petz2012efficient,petz2012optimal,rastegin2014notes,renes2004symmetric}.

    Our proposal is based on the application of the Naimark's theorem. The POVM can't be realized by an unitary transformation followed by a projective measurement in the original qudit space. Therefore, the initial qudit Hilbert space is expanded and the realization of the POVM elements occurs by applying unitary transformations followed by projection measurements in this extended space  \cite{Barnett,beneduci2010infinite,coecke2008povms,beneduci2010unsharpness,tabia2012experimental}. This task can be realized in multiport interferometers, composed of beam splitters and phase shifters, built by using, for example, the Reck's proposition  \cite{Reck94}. An important advance was done by Clements \textit{et al.}, that propose a different design for the multiport interferometer that, by using the same number of components, decreases the optical depth of the interferometer \cite{clements2016optimal}. This quantity is important for fabricating smaller interferometers with smaller losses.

    The unitary transformations in the extended Hilbert space in our circuit design uses a smaller number of beam splitters than Reck \textit{et al.} and Clements \textit{et al.} general designs with a smaller optical depth, what is a clear advantage for this specific task. It is necessary to mention that the proposed circuit design realizes full quantum tomography with the minimum number of measurements and no hypothesis about the initial system state are necessary. Others schemes are able to realize quantum tomography with a number of measurements smaller than the minimum for a particular class of systems states using prior knowledge of the input state \cite{oren2017quantum}.

    Finally, the organization of the rest of the paper is as follows: in Sec.~\ref{sec:tomography}, we briefly discuss some useful concepts about minimum state tomography for qubits. After this short introduction, we present in Sec.~\ref{sec:circuitsqubits} the proposal for construction of the tomographical circuits for this case. We then proceed in Sec.~\ref{sec:circuitsqutrit} with the circuit's proposal in the qutrit case. In Sec.~\ref{sec:circuitsNd}, we generalize our proposal for $N$ dimensions. Finally, we conclude in Sec.~\ref{sec:conclusao}.

	\section{Minimum quantum state tomography for qubits}
	\label{sec:tomography}

    As already mentioned in the introduction, quantum state tomography is an experimental technique that allows the reconstruction of the density matrix of an unknown quantum system state. Its realization consists in the production of a large number of identically prepared system states together with a series of measurements of the physical quantities that it describes. When the tomography is made using the least possible number of operations, it is called a minimum full quantum state tomography \cite{Rehacek, Pimenta10}. To perform a minimal full tomography on a qubit system, one needs a POVM with four different elements $\hat{E}_{j}$ to be implemented in the initial state $\hat{\rho}$. The outcome's counting rate $c_j$ is
     \begin{equation}
	    c_{j}=AP_j,
	    \label{count}
	\end{equation}
	where $A$ is a constant proportional to the detector efficiency and the number of identically prepared states in one second and $P_j=\text{Tr}[\hat{\rho}\hat{E}_{j}]$ is the probability of detecting one of the output states after $E_j$ is implemented. By implementing these POVM elements, it is possible to determine the qubit state since the number of counts is going to be accessible.

	As seen in \cite{Pimenta10}, the necessary POVM elements to be implemented are:
	\begin{equation} \label{POVMs}
		\hat{E}_{k}=\dfrac{1}{2}\ket{\varphi_{k}}\bra{\varphi_{k}}
	\end{equation}
	with $k=1,2,3,4$ and
	\begin{subequations}
		\begin{align}
		\ket{\varphi_1} & = \omega \ket{0} + i \upsilon \ket{1}, \\
		\ket{\varphi_2} & = \omega \ket{0} - i \upsilon \ket{1}, \\
		\ket{\varphi_3} & = \upsilon \ket{0} - \omega \ket{1}, \\
		\ket{\varphi_4} & = \upsilon \ket{0} + \omega \ket{1},
		\end{align}
	\end{subequations}
	where $\omega=\sqrt{2/3}$ and $\upsilon=\sqrt{1/3}$. If the density operator $\hat{\rho}$ is written as $\hat{\rho}=(\hat{I}+\vec{r}\cdot\vec{\sigma})/2$, where $\vec{r}=(r_x, r_y, r_z)$ and $\vec{\sigma}=(\hat{\sigma}_x,\hat{\sigma}_y,\hat{\sigma}_z)$ is the Pauli vector, we need to determine the values of the quantities $r_x$, $r_y$ and $r_z$ from the experimental outcome's counting rates $c_j$ (Eq.~\eqref{count}) in order to reconstruct the initial state density matrix $\hat{\rho}$. By using the POVM elements shown in Eq.~\eqref{POVMs} in Eq.~\eqref{count}, we obtain
	\begin{subequations}
		\begin{align}
			r_x=\dfrac{3(c_4-c_3)}{A\sqrt{2}},\qquad \\
		r_y=\dfrac{3(c_1-c_2)}{A\sqrt{2}}, \qquad	\\
		r_z=\dfrac{3(c_1+c_2-c_3-c_4)}{A},
		\end{align}
	\end{subequations}
	where $A=c_1+c_2+c_3+c_4$.

	\section{Construction of the tomographical circuits for qubits}
	\label{sec:circuitsqubits}

    Any unitary operator of finite size can be constructed in the laboratory using an appropriate combination of two optical devices: beam splitters (BS) and phase shifters \cite{Reck94}. Therefore, we can express an $N$-dimensional operator $\hat{U}'$ as:
    \begin{equation}\label{circuitomatrix}
	   \hat{U}'= \hat{T}_n \cdot \hat{T}_{n-1} \cdots \hat{T}_2 \cdot \hat{T}_1,
	\end{equation}
    where $\hat{T}_j$ is an $N$-dimensional matrix that represents the $j$th operation performed on the input state and $n$ is the number of operations required for the implementation of the matrix $\hat{U}'$. The $\hat{T}_j$ matrix is a block-diagonal matrix, where each block is composed from the elements of the BS matrix, defined as follows \cite{Reck94}:
    \begin{equation}
	    \hat{U}_{BS} =
	    \left[
	    \begin{array}{cc}
	        e^{i\phi}\sqrt{r}  & e^{i\phi}\sqrt{t} \\
	        \sqrt{t}           & -\sqrt{r}
	    \end{array}
	    \right],
	\end{equation}
    in which $t$ is the transmittance of the beam splitter and $r$ is its reflectivity. These parameters satisfy the equation $r+t=1$. The parameter $\phi$ indicates the phase difference between the outputs of the beam splitter. These are the building blocks we will use to create a scheme for minimal quantum state tomography in integrated chips.

    The first goal of this article is to propose a photonic circuit that implements the minimal full path state tomography on an one-qubit system. The simpler choice for the interferometer is shown in Fig.~\ref{saidaspovm}. In this configuration, each outcome of the interferometer corresponds to the implementation of one of the different maps $\hat{E}_k$, with $k=1,2,3,4$ (Eq.~\eqref{POVMs}), allowing us to obtain the outcomes simultaneously. The $\hat{T}_j$ operators that it constitutes (Fig.~\ref{saidaspovm}) are:
	\begin{equation}
	    \hat{T}_1 =
	    \left[
	    \begin{array}{cccc}
	        e^{i \phi_{1}} \sqrt{r_1} & e^{i \phi_{1}} \sqrt{t_{1}} & 0                          & 0 \\
	        \sqrt{t_{1}}                & -\sqrt{r_{1}}               & 0                          & 0 \\
	        0                           & 0                           & e^{i\phi_{2}} \sqrt{r_{2}} & e^{i\phi_{2}} \sqrt{t_{2}} \\
	        0                           & 0                           & \sqrt{t_{2}}               & -\sqrt{r_{2}}
	    \end{array}
	    \right],
	\end{equation}
	\begin{equation}
	    \hat{T}_2 =
	    \left[
	    \begin{array}{cccc}
	        1 & 0                           & 0                           & 0 \\
	        0 & e^{i \phi_{3}} \sqrt{r_{3}} & e^{i \phi_{3}} \sqrt{t_{3}} & 0 \\
	        0 & \sqrt{t_{3}}                & -\sqrt{r_{3}}               & 0 \\
	        0 & 0                           & 0                           & 1
	    \end{array}
	    \right],
	\end{equation}
	\begin{equation}
	\hat{T}_3 =
	\left[
	\begin{array}{cccc}
	e^{i \phi_{4}} \sqrt{r_4} & e^{i \phi_{4}} \sqrt{t_{4}} & 0                          & 0 \\
	\sqrt{t_{4}}                & -\sqrt{r_{4}}               & 0                          & 0 \\
	0                           & 0                           & e^{i\phi_{5}} \sqrt{r_{5}} & e^{i\phi_{5}} \sqrt{t_{5}} \\
	0                           & 0                           & \sqrt{t_{5}}               & -\sqrt{r_{5}}
	\end{array}
	\right].
	\end{equation}

    Since we need to expand the Hilbert space, we will use as initial state $\ket{\psi}_0=\ket{\psi}\otimes\ket{0}$. The first state vector can be interpreted as a qubit system in the state $\ket{\psi}$, and the second as an \textit{ancilla} in the state $\ket{0}$. But, note that the state $\ket{\psi}_0=(\alpha,0,\beta,0)^T$ is different from the initial state shown in the Fig.~\ref{saidaspovm}, given by $\ket{\psi}_{in}=(0,\alpha,\beta,0)^T$. The input state $\ket{\psi}_{in}$ can be written in terms of the state $\ket{\psi}_0$, with no modifications in the circuit nor any relabeling, by using the transformation described by $\hat{M}$, defined by the following relation
    \begin{equation}
    \left[
    \begin{array}{cccc}
    0 \\ \alpha \\ \beta \\  0
    \end{array}
    \right]=\hat{M}\cdot
    \left[
    \begin{array}{cccc}
    \alpha \\ 0 \\ \beta \\  0
    \end{array}
    \right]
    \quad\Rightarrow\quad
    \label{M}
	    \hat{M} =
	    \left[
	    \begin{array}{cccc}
		    0 & 1 & 0 & 0 \\
		    1 & 0 & 0 & 0 \\
		    0 & 0 & 1 & 0 \\
		    0 & 0 & 0 & 1
	    \end{array}
	    \right].
    \end{equation}

    \begin{figure}[]
    	\centering
    	\hspace*{-0.2cm} \includegraphics[height=4cm]{Fig_1}
    	\caption{Circuit proposal for the realization of the minimal state tomography on an one-qubit system. The black lines are the photon paths realized by the waveguides; the BS's are represented by the approximation of these lines and the phase shifters are the blue rectangles. This proposal was made on the assumption that the modes of the photonic qubit labelled as $\ket{0}$ and $\ket{1}$ are coupled to the second and third waveguides, respectively.}
    	\label{saidaspovm}
    \end{figure}

    It is worth noting that the transformation $\hat{M}$ is not implemented in practice. It is only a mathematical transformation. Thus, as our proposal supposes the input state as being  $\ket{\psi}_{in}$ and the theoretical calculation is made based on the initial state $\ket{\psi}_0$, we have
   \begin{equation}
   \hat{U}'\ket{\psi}_{in} = \hat{T}_3 \cdot \hat{T}_2 \cdot \hat{T}_1 \cdot \hat{M}\ket{\psi}_0.
   \end{equation}

   Therefore, the operator that effectively would act on the state $\ket{\psi}_0$ is given by
    \begin{equation}\label{portaU}
	    \hat{U}_{qubit} = \hat{T}_3\cdot\hat{T}_2\cdot\hat{T}_1\cdot\hat{M}.
	\end{equation}

	The operator defined in Eq.~\eqref{portaU} is used in the process of determination of the circuit parameters, as shown in the Appendix A.

	The next task is to find the parameters $r_{j}$, $t_{j}$ and $\phi_{j}$ that compose $\hat{U}_{qubit}$ to complete the circuit description.
	We start by examining the set of probability operators $\hat{\Pi}_{ml}$ acting on the system-ancilla Hilbert space. Such operators can be understood as the POVM elements associated with the detection of the photon at each circuit output. The action results of this POVM elements correspond to the projections over the states of the set $\{\ket{ml}\}$, with $m,l=0,1$, and each circuit output corresponds to one of the states $\ket{ml}$, i.e., $\ket{00}$, $\ket{01}$, $\ket{10}$ or $\ket{11}$. Here, the states $\ket{m}$ and $\ket{l}$ form a basis on each of the subspaces: system and \textit{ancilla}, respectively. The probability of detecting a photon, initially in the state described by the density operator $\hat{\rho}$, in one given output of the interferometer is given by
	\begin{equation}\label{probml}
		P_{ml}=\text{Tr}\left[\left(\hat{\rho}\otimes\ket{0}\bra{0}\right)\hat{\Pi}_{ml}\right],
	\end{equation}
    where $\ket{0}$ is the \textit{ancilla} state. This expression can also be written in terms of a POVM acting only in the subspace of the input qubit
	\begin{equation}
	P_{ml} = \text{Tr} \left[ \hat{\rho} \hat{E}_{ml} \right],
	\end{equation}
	where we can define the local POVM elements $\hat{E}_{ml}$ as
	\begin{equation} \label{traco2}
	\hat{E}_{ml} = \text{Tr}_2 \left[ \left( \hat{\mathbb{I}} \otimes \ket{0}\bra{0} \right) \hat{U}_{qubit}^\dagger\ket{ml}\bra{ml} \hat{U}_{qubit} \right]
	\end{equation}
    and the matrix elements of $\hat{E}_{ml}$ will be given by
    \begin{equation} \label{barneteq}
    \bra{p}\hat{E}_{ml}\ket{q} = \bra{p0} \hat{U}_{qubit}^{\dagger} \ket{ml} \bra{ml} \hat{U}_{qubit} \ket{q0}.
    \end{equation}

	Each state $\ket{ml}$, in binary notation, corresponds to one outcome shown in Fig.~\ref{saidaspovm}, which allows us to make the exchange $\hat{E}_{ml}\rightarrow\hat{E}_k$, being $k$ the number that identifies each output. Based on this, the operators $\hat{E}_{k}$ in Eq.~\eqref{traco2} will be identified as the POVM element presented in Eq.~\eqref{POVMs}. From these relations, it is possible to determine the suitable circuit parameters.
	As shown in the Appendix A, the parameters necessary for implementing the minimal tomography, in this case, are:
	\begin{subequations}
		\begin{align}
		r_1 & = 1/3, & t_1 & = 2/3, & \phi_1 & = \pi/2, \\
		r_2 & = 1/3, & t_2 & = 2/3,  & \phi_2 & = 0, \\
		\label{swap}
		r_3 & = 0,  & t_3 & = 1, & \phi_3 & = 0, \\
		r_4 & = 1/2, & t_4 & = 1/2, & \phi_4 & = 0, \\
		r_5 & = 1/2, & t_5 & = 1/2, & \phi_5 & = 0.
		\end{align}
	\end{subequations}

	By manufacturing a photonic circuit (Fig.~\ref{saidaspovm}) with the parameters presented in Eq. (17), we obtain a specific interferometer with five beam splitters and one phase shifter. This circuit implements POVM by doing unitary transformations and projective measurements in a extended Hilbert space with less optical elements and a smaller optical depth than other existing general unitary operations circuits already proposed in the literature \cite{Reck94,clements2016optimal,tabia2012experimental}. Other conclusion that arises from this result, is that the effect of the beam splitter that composes $\hat{T}_2$ solely exchanges the photon beam carried by the 2nd and 3rd paths, as seen in Eq.~\eqref{swap}. Except by this exchange, the BS's that acts non-trivially operates only between the 1st and 2nd, and the 3rd and 4th, separately. This structure will inspire the design of circuits for general qudit states.

	\section{Minimum quantum state tomography for qutrits}
	\label{sec:circuitsqutrit}

	For the qutrit, we approached the circuit's design in a different manner. The strategy to find a circuit that perform qubit's tomography starts with Eq.~\eqref{portaU}, which lead to the solution of a $4\times4$ matrix equation for finding the BS specifications. Applying the same procedure in the qutrit case would result in the solution of a $9 \times 9$ matrix equation. Alternatively, we faced this problem by analyzing the POVM operators and the desired outcomes for the qutrit tomography, being able to reduce the problem to the search of a $3 \times 3$ matrix equation solution. The qutrit tomography is realized via equidistant states \cite{Pimenta13,PaivaSanchez}.

	The SIC-POVM elements for the qutrit case are of the form
	\begin{equation}\label{purestates}
	    \hat{E}_{ml} = \frac{1}{3}| \psi_{ml}\rangle \langle \psi_{ml} | ,
	\end{equation}
	with $m,l=0,1,2$ and \cite{PaivaSanchez}
	\begin{subequations}
		\begin{align}
		|\psi_{00}\rangle & = \mu| 0 \rangle + \eta | 2 \rangle,\\
		|\psi_{01}\rangle & = \kappa| 0 \rangle + \kappa | 2 \rangle,  \\
		|\psi_{02}\rangle & = \eta| 0 \rangle + \mu | 2 \rangle,  \\
		|\psi_{10}\rangle & = \eta| 0 \rangle + \mu | 1 \rangle,   \\
		|\psi_{11}\rangle & = \kappa| 0 \rangle + \kappa | 1 \rangle, \\
		|\psi_{12}\rangle & = \mu| 0 \rangle + \eta | 1 \rangle,  \\
		|\psi_{20}\rangle & = \eta| 1 \rangle + \mu | 2 \rangle,  \\
		|\psi_{21}\rangle & = \kappa| 1 \rangle + \kappa | 2 \rangle,  \\
		|\psi_{22}\rangle & = \mu| 1 \rangle + \eta | 2 \rangle,
		\end{align}
	\end{subequations}
	 where $\kappa = 1/\sqrt{2}$, $\mu= e^{2\pi i / 3}/\sqrt{2}$, and $\eta = e^{- 2\pi i / 3}/\sqrt{2}$, such that $\kappa + \mu + \eta = 0$, ensuring that $\sum_{ml} \hat{E}_{ml} = \mathds{1}$. Keeping in mind Eq.~\eqref{barneteq}, we use the fact that the POVM elements are of the form of Eq.~\eqref{purestates} to define $a^{p}_{ml}$ as the expansion coefficients of $\ket{\psi_{ml}}$ as
	\begin{equation} \label{psiml}
	    \frac{1}{\sqrt{3}}\ket{\psi_{ml}} = a^{0}_{ml} \ket{0} + a_{ml}^1 \ket{1} + a_{ml}^2 \ket{2},
	\end{equation}
	where $\bra{p} \hat{E}_{ml} \ket{q} = (a_{ml}^p)^* a_{ml}^q$, with $p,q=0,1,2$. By substituing Eq.~\eqref{psiml} in Eq.~\eqref{purestates} and comparing $\hat{E}_{ml}$ with Eq.~\eqref{barneteq}, we arrive at an expression for the matrix elements of $\hat{U}$ that allow us to implement the desired POVM elements in the extended Hilbert space of dimension 9
	\begin{equation}
	\bra{ml} \hat{U} \ket{q0} = a_{ml}^q.
	\end{equation}

	\begin{figure}[]
		\centering
		\includegraphics[height=5.7cm]{Fig_2}
		\caption{Schematic drawing showing the quantum operations of the proposed photonic circuit for quantum tomography of qutrits. The pink dotted lines separate the three sectors of the circuit, that can be analyzed separately. The waveguides are represented by the black lines and labeled by the vectors $\ket{i}_j$, where $i,j=0,1,2$. The state vectors $\alpha_i\ket{i}_1$ indicate the waveguides where a possible $(\alpha_0,\alpha_1,\alpha_2)^T$ input qutrit path state is coupled to the circuit. The $\hat{\Lambda}$ operation is responsible for diffusing the initial state to the first and third path in each sector. The red and blue arrows portray the permutations between the $3$rd and $4$th paths and between the $6$th and $7$th ones, respectively. This permutation, performed by $1:0$ BS, exchanges the coefficients of some base vectors. Finally, these paths are connected to the last circuit piece, described by a $3\times3$ $\hat{U}$ matrix, to resume the tomographic implementation. Analogous descriptions can be done for the others two sectors. A detector count in one of the nine circuit outputs is proportional to the implementation probability of a POVM element $\hat{E}_{ij}$ $(i,j=0,1,2)$.}
		\label{pedacoqutrit}
	\end{figure}

	\begin{figure}[]
		\centering
		\includegraphics[height=4cm]{Fig_3}
		\caption{Circuit representation of the $\hat{\Lambda}$ operation for the implementation of the qutrit tomography. We have a $1/2:1/2$ BS  that splits the input mode in path 2 in two output modes in paths 2 and 3 (labelled as $\ket{i}_1$ and $\ket{i}_2$, respectively). A $1:0$ BS between paths 1 and 2 switches these paths.}
		\label{qutrit}
	\end{figure}

	Inspired by the solution for the tomography of qubits, we conceive an interferometer that is formed by three sectors, shown in Fig.~\ref{pedacoqutrit}, separated by the pink dotted line. Due to the symmetry in the POVM, it is convenient to perform the tomography measurements, such that, in the $k$th sector, the three outcomes result from the application of only $\hat{E}_{kl}$ $(l=0,1,2)$, already defined in Eq.~\eqref{purestates}. Our goal is that the photon count in the detector coupled to one of the three exits of the zeroth, first and second sector is proportional to the probability of implementing the POVM elements $\hat{E}_{0l}$, $\hat{E}_{1l}$ or $\hat{E}_{2l}$, respectively.

	The prepared photonic path state is such that the photon can enter in one of the nine input ports of the interferometer, as shown in Fig.~\ref{pedacoqutrit}. A photon in a qutrit path state represented by the $\ket{0}$, $\ket{1}$ and $\ket{2}$ base vectors can enter in the second or in the fifth or in the eighth input circuit ports. Since the photonic circuit has nine input ports, it is clear that the Hilbert space for the photonic path states increases from $N=3$ to $N=9$. We label the path states in the $N=9$ extended Hilbert space as $\ket{i}_j$, where $i,j=0,1,2$. The nine photonic paths, from top to bottom, are labelled as: $\ket{0}_0$,$\ket{0}_1$,$\ket{0}_2$, $\ket{1}_0$,$\ket{1}_1$,$\ket{1}_2$, $\ket{2}_0$,$\ket{2}_1$ and $\ket{2}_2$. Fig.~\ref{pedacoqutrit} shows the quantum operations schematically in the extended Hilbert space.

	Let's suppose we have a normalized input state $\ket{\psi}=\sum_{i=0}^{2}\alpha_i\ket{i}_1$. After the first operation, represented by the green rectangles, the photon path state becomes $\ket{\psi'}=\sum_{i=0}^{2}\alpha_i\left(\ket{i}_0+\ket{i}_2\right)/\sqrt{2}$. The second transformation in Fig.~\ref{pedacoqutrit} switches some of the coefficients in $\ket{\psi'}$ to generate the state $\ket{\psi''}=\sum_{i=0}^{2}\left(\beta_i\ket{i}_0+\gamma_i\ket{i}_2\right)/\sqrt{2}$, where $\beta_0=\beta_1=\alpha_0$, $\beta_2=\gamma_0=\alpha_1$ and $\gamma_1=\gamma_2=\alpha_2$. The states $\ket{\psi}$ and $\ket{\psi''}$ can be represented by vectors such that the first and second state transformations in Fig.~\ref{pedacoqutrit} produces

	\begin{equation}
		\left[
		\begin{array}{ccc}
		0 \\ \alpha_0 \\ 0 \\ 0 \\ \alpha_1 \\ 0 \\ 0 \\ \alpha_2 \\ 0 \\
		\end{array}
		\right] \quad\longmapsto\quad \dfrac{1}{\sqrt{2}}
		\left[
		\begin{array}{ccc}
		\alpha_0 \\ 0 \\ \alpha_1 \\ \alpha_0 \\ 0 \\ \alpha_2 \\ \alpha_1 \\ 0 \\ \alpha_2 \\
		\end{array}
		\right].
	\end{equation}

	\begin{figure}[]
	\centering
	\includegraphics[height=2.1cm]{Fig_4}
	\caption{Circuit representation of the $\hat{U}$ matrix for the implementation of the qutrit tomography. The blue rectangles are phase shifters.}
	\label{figuraU}
\end{figure}

	Notice that the circuit design is such that, in the first, second, and third sectors, we create state superpositions of the original qutrit base states $\{\ket{0}$, $\ket{1}\}$, $\{\ket{0}$, $\ket{2}\}$ and $\{\ket{1}$, $\ket{2}\}$, respectively. The arrangement of beam splitters which implements what is outlined by the green rectangles in Fig.~\ref{pedacoqutrit} is shown in Fig.~\ref{qutrit}. The optical couplers (BS) that carry the quantum operation described in Fig.~\ref{qutrit} are $1/2:1/2$ BS ($t=0.5$ and $r=0.5$) and $1:0$ BS ($t=1$ and $r=0$). The second transformation shown in Fig.~\ref{pedacoqutrit}, the permutations between 3rd and 4th paths and 6th and 7th paths, is performed by $1:0$ BSs.

	The last part of Fig.~\ref{pedacoqutrit} shows the applications of unitary operations $\hat{U}$ in the first, second and third sectors. Notice that the $\hat{U}$ operations only apply to the vector states of their respective sectors. Therefore, the quantum operation described by a $9\times9$ matrix in the extended Hilbert space of 9 dimensions is composed by three $3\times3$ $\hat{U}$ matrices forming three blocks matrices with all the other elements of the $9\times9$ matrix equal to zero, i.e.,
	\begin{equation}
		\hat{O}=\left[
		\begin{array}{ccc}
		\hat{U} &  0 & 0 \\
		 0 & \hat{U} & 0 \\
		 0 & 0 & \hat{U}  \\
		\end{array}
		\right].
	\end{equation}

The circuit scheme for $\hat{U}$ appears in detail in Fig.~\ref{figuraU}. Hence, we must establish the coefficients $r_j$, $t_j$, and the phases $\phi_j$. In order to achieve this, we are going to obtain the matrix from the inputs and outputs vectors available. The desired effect of the operations $\hat{E}_{0l}$, $\hat{E}_{1l}$ and $\hat{E}_{2l}$ in the first, second and third sectors of the circuit, respectively, leads us to
\begin{equation}
	\hat{U}\left[
	\begin{array}{ccc}
	\alpha_i/\sqrt{2} \\ 0 \\ \alpha_j/\sqrt{2}
	\end{array}
	\right]=\dfrac{1}{\sqrt{6}}
	\left[
	\begin{array}{ccc}
	\alpha_i+\alpha_j \\ \alpha_ie^{2i\pi/3}+\alpha_je^{-2i\pi/3} \\ \alpha_ie^{-2i\pi/3}+\alpha_je^{2i\pi/3}
	\end{array}
	\right],
	\label{u}
\end{equation}
with $i=0,1$, $j=1,2$ and $i<j$. The unitary $\hat{U}$ that satisfies Eq.~\eqref{u} is
	\begin{equation}
	\hat{U} = \frac{1}{\sqrt{3}}
	\left[
	\begin{array}{ccc}
	1       & 1       & 1 \\
	e^{2i\pi/3}   & 1   & e^{-2i\pi/3} \\
	e^{-2i\pi/3}       & 1       & e^{2i\pi/3}
	\end{array}
	\right],
	\label{operana}
	\end{equation}
	that corresponds to the circuit in Fig.~\ref{figuraU}. We managed therefore to reduce the problem from a $9 \times 9$ matrix equation to a problem of finding coefficients for a $3 \times 3$ matrices. Performing this calculation (more details in Appendix B), we obtain the values of the coefficients and phases:
	\begin{subequations}
		\begin{align}
		r_1 & = 1, & t_1 & = 0, & \phi_1 & = -2\pi/3, \\
		r_2 & = 0, & t_2 & = 1,  & \phi_2 & = -\pi/3, \\
		r_3 & = 1/2,  & t_3 & = 1/2, & \phi_3 & = -\pi/2, \\
		r_4 & = 1/3, & t_4 & = 2/3, & \phi_4 & = 0, \\
		r_5 & = 1/2, & t_5 & = 1/2, & \phi_5 & = \pi.
		\end{align}
	\end{subequations}

	The result presented in Eq. (26) shows that the first beam splitter of the scheme of Fig.~\ref{figuraU} can be discarded, since there is no transmittance. The phase shifter $\phi_5$ can also be discarded, since it is at the end of the circuit and does not interfere in the photon counting. We conclude that the interferometer that implements the quantum path state tomography in the one-qutrit state is formed by twenty beam splitters and nine phase shifters. As in the one-qubit case, our proposal for one-qutrit systems has less optical elements and a smaller optical depth than other existing proposals in the literature \cite{Reck94,clements2016optimal,tabia2012experimental}.

	We also summarize in Table~\ref{table} the number of optical elements necessary for the implementation of each proposal of photonic circuits for the quantum tomography task. It is worth emphasizing that Reck's and Clements's proposals use Mach-Zehnder interferometers as building blocks of the circuits, while our proposal uses single beam splitters.

	\begin{table}[]
		\caption{Number of optical elements required for the implementation of minimum quantum path state tomography in photonic circuits according to our proposal and others proposals already known in the literature\cite{Reck94,clements2016optimal,tabia2012experimental}.}
		\label{table}
		\begin{tabular}{|c|c|c|c|c|}
			\hline
			\multicolumn{1}{|l|}{} & \textbf{\begin{tabular}[c]{@{}c@{}}Reck's \\ proposal\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Clements's \\ proposal\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Tabia's \\ proposal\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Our \\ proposal\end{tabular}} \\ \hline
			\textbf{Qubit}         & 12                                                                  & 12                                                                      & 7                                                                    & 6                                                                \\ \hline
			\textbf{Qutrit}        & 72                                                                  & 72                                                                      & 44                                                                   & 29                                                               \\ \hline
		\end{tabular}
	\end{table}

	\section{Generalization for $N$-dimensions}
	\label{sec:circuitsNd}

	In the previous section, we use the results of Paiva-S\'anchez \cite{PaivaSanchez} about equidistant states to build the SIC-POVM elements that implement a quantum state tomography of path photonic qutrits. One remarkable feature of this approach is that it is not limited to the dimension of the input state. The relation for equidistant states suitable for the tomography of an $N$-dimensional quantum system is \cite{PaivaSanchez}:
	\begin{equation}
	\ket{\psi_{l}}=\dfrac{1}{\sqrt{N-1}}\sum_{\substack{k=0 \\ k\neq 1}}^{N-1}e^{2i\pi l\left(k-1\right)/N}\ket{k},
	\label{equidistante2}
	\end{equation}
	where $l=0,1,...,N-1$. From Eq.~\eqref{equidistante2}, we can to obtain $N$ operators of the $N^2$ required for the quantum state tomography experiment. The others operators are obtained from states generated by the application of the operator $\hat{X}$ in $\ket{\psi_l}$ shown in Eq.~\eqref{equidistante2}, defined as:
	\begin{equation}
	\hat{X}\ket{k}=\ket{k\oplus 1}.
	\label{x}
	\end{equation}
	This operator $\hat{X}$ acts on a state $\ket{k}$ performing an addition modulo $N$. Thus, all equidistant states used in the construction of the POVM elements can be organized in sets $B_m(\psi)$, with $m=0,1,...,N-1$, defined as
	\begin{equation}
	B_m(\psi)=\left\{\ket{\psi_{ml}}=\hat{X}^m\ket{\psi_l}\right\}
	\end{equation}
	and these POVM elements are obtained through the relation
	\begin{equation}
\hat{E}_{ml} = \dfrac{1}{N}| \psi_{ml}\rangle \langle \psi_{ml}|.
\end{equation}

		As an example, we will define the base states for quantum state tomography of ququarts, that is, qudits with $N=4$. They are:
	\begin{subequations}
	\begin{align}
		\ket{\psi_{00}}&=\dfrac{1}{\sqrt{3}}\left(\ket{0}+\ket{2}+\ket{3}\right); \\
		\ket{\psi_{01}}&=\dfrac{1}{\sqrt{3}}\left(-i\ket{0}+i\ket{2}-\ket{3}\right); \\
		\ket{\psi_{02}}&=\dfrac{1}{\sqrt{3}}\left(-\ket{0}-\ket{2}+\ket{3}\right); \\
		\ket{\psi_{03}}&=\dfrac{1}{\sqrt{3}}\left(i\ket{0}-i\ket{2}-\ket{3}\right).
	\end{align}
	\end{subequations}

	The other states ($\ket{\psi_{1l}}$, $\ket{\psi_{2l}}$ and $\ket{\psi_{3l}}$, with $l=0,1,2,3$) come from the application of the operator $\hat{X}$ one, two or three times.

	The determination of the circuit design for the general case is made in the same way as it was done in the Section IV, not applicable for qubit case. The $N^2$ paths are divided in the $N$ sectors with $N$ paths each. Here, the three-dimensional conception of circuits facilitates the generalization task \cite{chaboyer2015tunable,meany2015laser,poulios2014quantum}. The assumption of a 3D waveguide writing technique becomes necessary because we have to create path states superposition between states of different sectors. For example, it is necessary to create superposition between the states of the first and the last sector in a SIC-POVM tomographic proposal. In a planar circuit, the optical depth would increase, as long as these states are connected, by a quadratic function in $N$. In the 3D waveguide writing technique, these states will be connected by $N-1$ beam splitters, that is, the optical depth is now a linear function in $N$.

	In the three-dimensional circuit scheme, the paths can be organized in $N$ parallel vertical sectors, with $N$ paths in the vertical direction each, forming a square array of waveguides. By choosing the paths for entrance of the photons belonging to the same line, i.e., path $i_0$ in each sector $j$, we have as initial state
	\begin{equation}
		\ket{\psi}_j^{(0)}=\sum_{j=0}^{N-1}\alpha_{j}\ket{i_0j},
		\label{initial}
	\end{equation}
	where $\ket{i_0j}$ means path state in line $i_0$ and sector $j$. Our proposal consists in conceiving the circuit as a sequence of three unitary operations. The first one operates separately in each vertical sector $j$ performing the transformation
	\begin{equation}
		\alpha_{j}\ket{i_0j} \quad\longmapsto\quad \dfrac{\alpha_{j}}{\sqrt{N-1}}\sum_{i=0}^{N-2}\ket{ij},
	\end{equation}
	i.e., each component of the initial state $\alpha_j\ket{i_0j}$ is decomposed into $N-1$ components $\alpha_j\ket{ij}$, that will occupy $N-1$ paths of the vertical sector $j$, leaving, also by choice, the last path empty.

	The second operation realized by the circuit is a permutation betweens vertical sectors. Except the last layer, that did not receive one of the $N-1$ components at the end of the first operation, in all other layers the operation performed is an addiction modulo $N$ as follows
	\begin{equation}
		\hat{\varGamma}\ket{ij}=\ket{i}\ket{j\oplus i},
	\end{equation}
	where $i=0,1,...,N-2$ is the paths label. At the end of this second operation, we will have in each vertical sector $j$ the state
	\begin{equation}
		\ket{\psi}_j^{(2)}=\dfrac{1}{\sqrt{N-1}}\sum_{i=0}^{N-2}\alpha_{j-i}\ket{ij},
		\label{psi2}
	\end{equation}
 	where the operations with subindexes are carried out modulo $N$. The third operation that our proposed circuit performs is analogous to Eq.~\eqref{operana}. Its transforms the state $\ket{\psi}_j^{(2)}$ (Eq.~\eqref{psi2}) in a final state where the square module of its expansion coefficients are igual to the probabilities of implement all POVM elements in the qudit Hilbert space. By possession of these probabilities, we reconstruct the initial state density matrix \cite{PaivaSanchez}.

	If the possibility of usage of a 3-dimensional array is at hand, our circuit becomes considerably more efficient that one derived directly from \cite{Reck94} or \cite{clements2016optimal}. As the minimal state tomography uses at least $N^2$ outputs, a general circuit would have a complexity of $O(N^4)$ beam splitters. In our realization, the procedure is divide into three parts, each acting in parallel over at maximum $N$ paths. This guarantees a complexity of $O(N^3)$. Moreover, each part may simplify separately, since they act differently, allowing an increase of efficiency through the detalied analysis of each part. The modular property of the circuit is the main feature of our proposed implementation of the circuit for minimal state tomography.

	There is a particularity in the even-dimension case, reported in Ref. \cite{PaivaSanchez}. It is shown that the method exposed in this paper cannot determine the imaginary part of some elements of the system's density matrix $\hat{\rho}$. For that, one more step is needed, the application of the following operator:
	\begin{equation}
	\hat{\Xi}=\sum_{k=0}^{N/2-1}\ket{k}\bra{k}+i\sum_{N/2}^{N-1}\ket{k}\bra{k},
	\label{ureal}
	\end{equation}
	which generates a new density matrix $\hat{\rho}'=\hat{\Xi}\hat{\rho}\hat{\Xi}^\dag$. The advantage about this strategy is that the imaginary parts, impossible to be determined, are exchanged with the real parts of the same matrix elements, being now possible to be obtained. Thus, to reconstruct the initial density matrix of even-dimensional systems, it is necessary to realize the quantum state tomography in $\hat{\rho}$ and $\hat{\rho}'$, increasing from $N^2$ to $3N^2/2$ the number of measurements.

	\section{Conclusion}
	\label{sec:conclusao}

	In this paper, we were able to determine the design of integrated photonic circuits that perform quantum path state tomography for qubits and qutrits. In these schemes, a POVM is implemented by doing an application of Naimark's theorem by extending the qudit Hilbet space, performing unitary transformations and projective measurements in this extended space. The parameters of the optical interferometer elements (beam splitters and phase shifters) were accurately calculated by using methods already described in the literature. Only three different types of beam splitters appeared on the proposed circuits, which may facilitate their manufacture. Our designs requires less optical components than alternatives general unitary transformations circuits and a smaller optical depth when compared with the optical depth of the circuits discussed in \cite{Reck94,clements2016optimal}. A general recipe also was presented for determining the circuit to be used for dimensions higher than three. In all cases, our interferometers are more compact than those already exposed in other works.

	\begin{acknowledgements} This research was supported by the Brazilian agencies CNPq - Conselho Nacional de Desenvolvimento Cient\'{\i}fico e Tecnol\'ogico, Capes - Coordena\c{c}\~ao de Aperfei\c{c}oamento de Pessoal de N\'{\i}vel Superior, Fapemig - Funda\c{c}\~ao de Amparo \`a Pesquisa do Estado de Minas Gerais and INCT-IQ - Instituto Nacional de Ci\^encia e Tecnologia de Informa\c{c}\~ao Qu\^antica.
	\end{acknowledgements}