text:
\section{Introduction}
\vspace{-0.1cm}
\begin{figure}[!t]
	\begin{subfigure}{.16\textwidth}
		\centering
		\includegraphics[width=\linewidth]{figures/mnist_test_pca-eps-converted-to.pdf}
		\caption{Raw Data}

	\end{subfigure}%
	\begin{subfigure}{.16\textwidth}
		\centering
		\includegraphics[width=\linewidth]{figures/mnist_test_sda_emc_pca-eps-converted-to.pdf}
		\caption{NonJoint \textit{DEPICT}}

	\end{subfigure}
	\begin{subfigure}{.153\textwidth}
		\centering
		\includegraphics[trim=40mm 30mm 30mm 30mm, clip,  width=1\linewidth]{figures/mnist_test_emced_pca_2}
		\caption{Joint \textit{DEPICT}}

	\end{subfigure}
	\caption{Visualization to show the discriminative capability of embedding subspaces using \textit{MNIST-test} data. (a) The space of raw data.  (b) The embedding subspace of non-joint \textit{DEPICT} using standard stacked denoising autoencoder (SdA). (c) The embedding subspace of joint \textit{DEPICT} using our joint learning approach (MdA).}
	\label{fig_1}
\end{figure}

Clustering is one of the fundamental research topics in machine learning and computer vision research, and it has gained significant attention for discriminative representation of data points without any need for supervisory signals. The clustering problem has been extensively studied in various applications; however, the performance of standard clustering algorithms is adversely affected when dealing with high-dimensional data, and their time complexity dramatically increases when working with large-scale datasets. Tackling the curse of dimensionality, previous studies often initially project data into a low-dimensional manifold, and then cluster the embedded data in this new subspace  \cite{roth2003feature,tian2014learning,wang2016learning}. Handling large-scale datasets, there are also several studies which select only a subset of data points to accelerate the clustering process \cite{shinnou2008spectral,Dijun:ICDE10,Huang2015AAAI2}.

However, dealing with real-world image data, existing clustering algorithms suffer from different issues: 1) Using inflexible hand-crafted features, which do not depend on the input data distribution; 2) Using shallow and linear embedding functions, which are not able to capture the non-linear nature of data; 3) Non-joint embedding and clustering processes, which do not result in an optimal embedding subspace for clustering; 4) Complicated clustering algorithms that require tuning the hyper-parameters using labeled data, which is not feasible in real-world clustering tasks.

To address the mentioned challenging issues, we propose a new clustering algorithm, called deep embedded regularized  clustering (\textit{DEPICT}), which exploits the advantages of both discriminative clustering methods and deep embedding models. \textit{DEPICT} generally consists of two  main parts, a multinomial logistic regression (soft-max) layer stacked on top of a multi-layer convolutional autoencoder. The soft-max layer along with the encoder pathway can be considered as a discriminative clustering model, which is trained using the relative entropy (\textit{KL} divergence) minimization. We further add a regularization term based on a prior distribution for the frequency of cluster assignments. The regularization term penalizes unbalanced cluster assignments and prevents allocating clusters to outlier samples.

Although this deep clustering model is flexible enough to discriminate the complex real-world input data,  it can easily get stuck in non-optimal local minima during training and result in undesirable cluster assignments. In order to avoid overfitting the deep clustering model to spurious data correlations, we utilize the reconstruction loss function of autoencoder models as a data-dependent regularization term for training parameters.

In order to benefit from a joint learning framework for embedding and clustering, we introduce a unified objective function including our clustering and auxiliary reconstruction loss functions. We then employ an alternating approach to efficiently update the parameters and estimate the cluster assignments. It is worth mentioning that in the standard learning approach for training a multi-layer autoencoder, the encoder and decoder parameters are first pretrained layer-wise using the reconstruction loss, and the encoder parameters are then fine-tuned using the objective function of the main task \cite{vincent2010stacked}. However, it has been argued that the non-joint fine-tuning step may overwrite the encoder parameters entirely and consequently cancel out the benefit of the layer-wise pretraining step \cite{zhao2015stacked}. To avoid this problem and achieve optimal joint learning results, we simultaneously train all of the encoder and decoder layers together along with the soft-max layer. To do so, we sum up the squared error reconstruction loss functions between the decoder and their corresponding (clean) encoder layers and add them to the clustering loss function.

Figure \ref{fig_1} demonstrates the importance of our joint learning strategy by comparing different data representations of \textit{MNIST-test} data points \cite{lecun1998gradient} using principle component analysis (PCA) visualization. The first figure indicates the raw data representation; The second one shows the data points in  the embedding subspace of non-joint \textit{DEPICT}, in which the model is trained using the standard layer-wise stacked denoising autoencoder (SdA); The third one visualizes the data points in the embedding subspace of joint \textit{DEPICT}, in which the model is trained using our multi-layer denoising autoencoder learning approach (MdA). As shown, joint \textit{DEPICT} using MdA learning approach provides a significantly more discriminative embedding subspace compared to non-joint \textit{DEPICT} using standard SdA learning approach.

Moreover, experimental results show that \textit{DEPICT} achieves superior or competitive results compared to the state-of-the-art algorithms on the image benchmark datasets while having faster running times. In addition, we compared different learning strategies for \textit{DEPICT}, and confirm that our joint learning approach has the best results. It should also be noted that  \textit{DEPICT} does not require any hyper-parameter tuning using supervisory signals, and consequently is a better candidate for the real-world clustering tasks.
Thus, we summarize the advantages of \textit{DEPICT} as:
\begin{itemize}[noitemsep,topsep=2pt]
	\item Providing a discriminative non-linear embedding subspace via the deep convolutional autoencoder;
	\item Introducing an end-to-end joint learning approach, which  unifies the clustering and embedding tasks, and avoids layer-wise pretraining;
	\item Achieving superior or competitive clustering results on high-dimensional and large-scale datasets with no need for hyper-parameter tuning using labeled data.
\end{itemize}

\section{Related Works}\label{section:realted_works}
There is a large number of clustering algorithms in literature, which can be grouped into different perspectives, such as hierarchical \cite{heller2005bayesian,williams1999mcmc,zhang2012graph}, centroid-based \cite{lloyd1982least,bezdek1984fcm,Huang2014KDD2,bahmani2012scalable}, graph-based \cite{shi2000normalized,nie2016constrained,wang2016structured,Huang2016AAAI5}, sequential (temporal) \cite{keogh2001online,sargin2008analysis,sadoughi2015msp,zhou2013hierarchical,sadoughi2015retrieving}, regression model based \cite{Huang2016ICDM1,Huang2013ICML1}, and subspace clustering models \cite{agrawal1998automatic,kailing2004density,Huang2015ICCV,Huang2016IJCAI1}. In another sense, they are generally divided into two subcategories, generative and discriminative clustering algorithms. The generative algorithms like \textit{$K$-means} and Gaussian mixture model \cite{biernacki2000assessing} explicitly represent the clusters using geometric properties of the feature space, and model the categories via the statistical distributions of input data. Unlike the generative clustering algorithms, the discriminative methods directly identify the categories using their separating hyperplanes regardless of data distribution. Information theoretic \cite{li2004minimum,barber2005kernelized,krause2010discriminative}, max-margin \cite{zhao2008efficient,xu2004maximum}, and spectral graph \cite{ng2002spectral} algorithms are examples of discriminative clustering models. Generally it has been argued that the discriminative models often have better results compared to their generative counterparts, since they have fewer assumptions about the data distribution and directly separate the clusters, but their training can suffer from overfitting or getting stuck in undesirable local minima \cite{krause2010discriminative,ng2002spectral,raina2003classification}. Our \textit{DEPICT} algorithm is also a discriminative clustering model, but it benefits from the auxiliary reconstruction task of autoencoder to alleviate this issues in training of our discriminative clustering algorithm.

There are also several studies regarding the combination of clustering with feature embedding learning. Ye \emph{et al.} introduced a kernelized \textit{$K$-means} algorithm, denoted by \textit{DisKmeans}, where embedding to a lower dimensional subspace via linear discriminant analysis (\textit{LDA}) is jointly learned with \textit{$K$-means} cluster assignments \cite{ye2008discriminative}. \cite{Huang2014ECML1} proposed to a new method to simultaneously conduct both clustering and feature embedding/selection tasks to achieve better performance. But these models suffer from having shallow and linear embedding functions, which cannot represent the non-linearity of real-world data.

A joint learning framework for updating code books and estimating image clusters was proposed in \cite{xie2015integrating} while \textit{SIFT} features are used as input data. A deep structure, named \textit{TAGnet} was introduced in \cite{wang2016learning}, where two layers of sparse coding followed by a clustering algorithm are trained with an alternating learning approach. Similar work is presented in \cite{wang2015joint} that formulates a joint optimization framework for discriminative clustering and feature extraction using sparse coding. However, the inference complexity of sparse coding forces the model in \cite{wang2015joint} to reduce the dimension of input data with \textit{PCA} and the model in \cite{wang2016learning} to use an approximate solution. Hand-crafted features and dimension reduction techniques degrade the clustering performance by neglecting the distribution of input data.

Tian \emph{et al.} learned a non-linear embedding of the affinity graph using a stacked autoencoder, and then obtained the clusters in the embedding subspace via \textit{$K$-means} \cite{tian2014learning}.
Trigeorgis \emph{et al.} extended semi non-negative matrix factorization (\textit{semi-NMF}) to stacked multi-layer (deep) \textit{semi-NMF} to capture the abstract information in the top layer.  Afterwards, they run \textit{$K$-means} over the embedding subspace for cluster assignments \cite{trigeorgis2014deep}.
More recently, Xie \emph{et al.} employed denoising stacked autoencoder learning approach, and first pretrained the model layer-wise and then fine-tuned the encoder pathway stacked by a clustering algorithm using Kullback-Leibler divergence minimization \cite{xie2016unsupervised}. Unlike these models that require  layer-wise pretraining as well as non-joint embedding and clustering learning, \textit{DEPICT} utilizes an end-to-end optimization for training all network layers simultaneously using the unified clustering and reconstruction loss functions.

Yang \emph{et al.} introduced a new clustering model, named \textit{JULE}, based on a recurrent framework, where data is represented via a convolutional neural network and embedded data is iteratively clustered using an agglomerative clustering algorithm \cite{yang2016joint}. They derived a unified loss function consisting of the merging process for agglomerative clustering and updating the parameters of the deep representation. While \textit{JULE} achieved good results using the joint learning approach, it requires tuning of a large number of hyper-parameters, which is not practical in real-world clustering tasks. In contrast, our model does not need any supervisory signals for hyper-parameter tuning.

\section{Deep Embedded Regularized Clustering} \label{section:method}
In this section, we first introduce the clustering objective function and the corresponding optimization algorithm, which alternates between estimating the cluster assignments and updating model parameters. Afterwards, we show the architecture of \textit{DEPICT} and provide the joint learning framework to simultaneously train all network layers using the unified clustering and reconstruction loss functions.

\subsection{DEPICT Algorithm} \label{section:RLC}

Let's consider the clustering task of $N$ samples, $\mathbf{X} = [\mathbf{x}_1, ..., \mathbf{x}_n]$, into $K$ categories, where each sample $\mathbf{x}_i \in \mathbb{R}^{d_x}$.  Using the embedding function, $\varphi_{W}: X\rightarrow Z$, we are able to map raw samples into the embedding subspace $\mathbf{Z} = [\mathbf{z}_1, ..., \mathbf{z}_n]$, where each $\mathbf{z}_i \in \mathbb{R}^{d_z}$ has a much lower dimension compared to the input data (i.e. $d_z \ll d_x$). Given the embedded features, we use a multinomial logistic regression (soft-max) function $f_{\theta}: Z\rightarrow Y$ to predict the probabilistic cluster assignments as follows.
\begin{align} \label{eq:1}
&	p_{ik} = P(y_i=k | \mathbf{z}_i, \boldsymbol{\Theta})= \frac{exp( \boldsymbol{\theta}_k^T \mathbf{z}_i )}{\sum\limits_{k'=1}^K exp ( \boldsymbol{\theta}_{k'}^T \mathbf{z}_i )} \,,
\end{align}
where $\boldsymbol{\Theta} = [\boldsymbol{\theta}_1, ..., \boldsymbol{\theta}_k] \in \mathbb{R}^{d_z \times  K}$ are the soft-max function parameters, and $p_{ik}$ indicates the probability of the $i$-th sample belonging to the $k$-th cluster.

In order to define our clustering objective function, we employ an auxiliary target variable $\mathbf{Q}$ to refine the model predictions iteratively. To do so, we first use Kullback-Leibler (\textit{KL}) divergence to decrease the distance between the model prediction $\mathbf{P}$ and the target variable $\mathbf{Q}$.
\begin{align} \label{eq:2}
\mathcal{L} = KL (\mathbf{Q} \| \mathbf{P}) = \frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  \frac{q_{ik}}{p_{ik}}  \,,
\end{align}
In order to avoid degenerate solutions, which allocate most of the samples to a few clusters or assign a cluster to outlier samples, we aim to impose a regularization term to the target variable. To this end, we first define the empirical label distribution of target variables as:
\begin{align} \label{eq:3}
&	f_{k} = P(\mathbf{y}=k) = \frac{1}{N} \sum\limits_{i} q_{ik}   \,,
\end{align}
where $f_k$ can be considered as the soft frequency of cluster assignments in the target distribution. Using this empirical distribution, we are able to enforce our preference for having balanced assignments by adding the following \textit{KL} divergence to the loss function.
\begin{align} \label{eq:4}
\mathcal{L} &= KL (\mathbf{Q} \| \mathbf{P})  + KL (\mathbf{f} \| \mathbf{u})  \\
&= \Big[ \frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  \frac{q_{ik}}{p_{ik}} \Big] +  \Big[ \frac{1}{N} \sum\limits_{k=1}^K f_{k} \log  \frac{f_{k}}{u_{k}} \Big] \nonumber \\
&= \frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  \frac{q_{ik}}{p_{ik}} +  q_{ik} \log  \frac{f_{k}}{u_{k}} \,, \nonumber
\end{align}
where $\mathbf{u}$ is the uniform prior for the empirical label distribution. While the first term in the objective minimizes the distance between the target and model prediction distributions, the second term balances the frequency of clusters in the target variables.
Utilizing the balanced target variables, we can force the model to have more balanced predictions (cluster assignments) $\mathbf{P}$ indirectly. It is also simple to change the prior from the uniform distribution to any arbitrary distribution in the objective function if there is any extra knowledge about the frequency of clusters.

An alternating learning approach is utilized to optimize the objective function. Using this approach,  we estimate the target variables $\mathbf{Q}$ via fixed parameters (expectation step), and update the parameters while the target variables $\mathbf{Q}$  are assumed to be known (maximization step). The problem to infer the target variable $\mathbf{Q}$ has the following objective:
\begin{align} \label{eq:5}
\min\limits_{\mathbf{Q}} \ \  \frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  \frac{q_{ik}}{p_{ik}} +  q_{ik} \log  \frac{f_{k}}{u_{k}} \,,
\end{align}
where the target variables are constrained to $\sum_k q_{ik} =1$. This problem can be solved using first order methods, such as gradient descent, projected gradient descent, and Nesterov optimal method \cite{nesterov2013introductory}, which only  require the objective function value and its (sub)gradient at each iteration. In the following equation, we show the partial derivative of the objective function with respect to the target variables.
\begin{align} \label{eq:6}
\frac{\partial  \mathcal{L}}{\partial q_{ik}} \propto  \log \Big( \frac{q_{ik} f_k}{p_{ik}} \Big) +  \frac{q_{ik}}{\sum\limits_{i'=1}^N q_{i'k}} + 1 \,,
\end{align}
Investigating this problem more carefully, we approximate the gradient in Eq.(\ref{eq:6}) by removing the second term, since the number of samples N is often big enough to ignore the second term. Setting the gradient equal to zero, we are now able to compute the closed form solution for $\mathbf{Q}$ accordingly.
\begin{align} \label{eq:7}
&	 q_{ik} =   \frac{  p_{ik}   / (\sum_{i'} p_{i'k})^{\frac{1}{2}}  }{\sum\limits_{k'}   p_{ik'}   / (\sum_{i'} p_{i'k'})^{\frac{1}{2}}}   \,,
\end{align}
For the maximization step, we update the network parameters $\boldsymbol{\psi} = \{\boldsymbol{\Theta}, \mathbf{W}\}$ using the estimated target variables with the following objective function.
\begin{align} \label{eq:8}
\min\limits_{\boldsymbol{\psi}} \ \  -\frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  p_{ik}\,,
\end{align}
Interestingly, this problem can be considered as a standard cross entropy loss function for classification tasks, and the parameters of soft-max layer $\boldsymbol{\Theta}$ and embedding function $\mathbf{W}$ can be efficiently updated by backpropagating the error.

\begin{figure*}[t]
	\begin{center}
		\includegraphics[trim=0mm 0mm 0mm 0mm, clip,  width=0.83\linewidth]{figures/drawing_last_no_rlc}
	\end{center}
	\vspace{-0.3cm}
	\caption{Architecture of \textit{DEPICT} for \textit{CMU-PIE} dataset. \textit{DEPICT} consists of a soft-max layer stacked on top of a multi-layer convolutional autoencoder. In order to illustrate the joint learning framework, we consider the following four pathways for \textit{DEPICT}: Noisy (corrupted) encoder, Decoder, Clean encoder and Soft-max layer. The clustering loss function, $L_E$, is applied on the noisy pathway, and the reconstruction loss functions, $L_2$, are between the decoder and clean encoder layers. The output size of convolutional layers,  kernel sizes, strides (S), paddings (P) and crops (C) are also shown. }
	\vspace{-0.2cm}
	\label{fig:2}
\end{figure*}

\subsection{DEPICT Architecture}
In this section, we extend our general clustering loss function using a denoising autoencoder. The deep embedding function is useful for capturing the non-linear nature of input data; However, it may overfit to spurious data correlations and get stuck in undesirable local minima during training. To avoid this overfitting, we employ autoencoder structures and use the reconstruction loss function as a data-dependent regularization for training the parameters. Therefore, we design \textit{DEPICT} to consist of a soft-max layer stacked on top of a multi-layer convolutional autoencoder. Due to the promising performance of strided convolutional layers in \cite{radford2015unsupervised,yeh2016semantic}, we employ convolutional layers in our encoder and strided convolutional layers in the decoder pathways, and avoid deterministic spatial pooling layers (like max-pooling). Strided convolutional layers allow the network to learn its own spatial upsampling, providing a better generation capability.

Unlike the standard learning approach for denoising autoencoders, which contains layer-wise pretraining and then fine-tuning, we simultaneously learn all of the autoencoder and soft-max layers. As shown in Figure \ref{fig:2}, \textit{DEPICT} consists of the following components:

\noindent \textbf{1)} Corrupted feedforward (encoder) pathway maps the noisy input data into the embedding subspace using a few convolutional layers followed by a fully connected layer. The following equation indicates the output of each layer in the noisy encoder pathway.
\begin{align} \label{eq:50}
&	 \tilde{\mathbf{z}}^{l} = Dropout \big[g(\mathbf{W}_e^l \tilde{\mathbf{z}}^{l-1} ) \big] \,,
\end{align}
where $\tilde{\mathbf{z}}^{l}$ are the noisy features of the $l$-th layer, $Dropout$ is a stochastic mask function that randomly sets a subset of its inputs to zero \cite{srivastava2014dropout}, $g$ is the activation function of convolutional or fully connected layers, and $\mathbf{W}_e^l$ indicates the weights of the $l$-th layer in the encoder. Note that the first layer features, $\tilde{\mathbf{z}}^0$, are equal to the noisy input data, $\tilde{\mathbf{x}}$.

\noindent \textbf{2)} Followed by the corrupted encoder, the decoder pathway reconstructs the input data through a fully connected and multiple strided convolutional layers as follows,
\begin{align} \label{eq:51}
&	 \hat{\mathbf{z}}^{l-1} = g(\mathbf{W}_d^l \hat{\mathbf{z}}^{l} )  \,,
\end{align}
where $\hat{\mathbf{z}}^{l}$ is the $l$-th reconstruction layer output, and $\mathbf{W}_d^l$ shows the weights for the $l$-th layer of the decoder. Note that input reconstruction, $\hat{\mathbf{x}}$, is equal to $\hat{\mathbf{z}}^0$.

\noindent \textbf{3)} Clean feedforward (encoder) pathway shares its weights with the corrupted encoder, and infers the clean embedded features. The following equation shows the outputs of the clean encoder, which are used in the reconstruction loss functions and obtaining the final cluster assignments.
\begin{align} \label{eq:52}
&	 \mathbf{z}^{l} = g(\mathbf{W}_e^l \mathbf{z}^{l-1} )  \,,
\end{align}
where $\mathbf{z}^{l}$ is the clean output of the $l$-th layer in the encoder. Consider the first layer features $\mathbf{z}^0$ equal to input data $\mathbf{x}$.

\noindent \textbf{4)} Given the top layer of the corrupted and clean encoder pathways as the embedding subspace, the soft-max layer obtains the cluster assignments  using Eq.(\ref{eq:1}).

Note that we compute target variables $\mathbf{Q}$ using the clean pathway, and model prediction $\tilde{\mathbf{P}}$ via the corrupted pathway. Hence, the clustering loss function $KL (\mathbf{Q} \| \tilde{\mathbf{P}})$ forces the model to have invariant features with respect to noise. In other words, the model is assumed to have a dual role: a clean model, which is used to compute the more accurate target variables; and a noisy model, which is trained to achieve noise-invariant predictions.

As a crucial point, \textit{DEPICT} algorithm provides a joint learning framework that optimizes the soft-max and autoencoder parameters together.
\begin{align} \label{eq:53}
\min\limits_{\boldsymbol{\psi}} \ \  -\frac{1}{N} \sum\limits_{i=1}^N \sum\limits_{k=1}^K q_{ik} \log  \tilde{p}_{ik}  + \frac{1}{N}  \sum\limits_{i=1}^N  \sum\limits_{l=0}^{L-1} \frac{1}{|\mathbf{z}_i^l|} \| \mathbf{z}_i^l - \hat{\mathbf{z}}_i^l \|_2^2 \,,
\end{align}
where $|\mathbf{z}_i^l|$ is the output size of the $l$-th hidden layer (input for $l=0$), and $L$ is the depth of the autoencoder model.

The benefit of joint learning frameworks for training multi-layer autoencoders is also reported in semi-supervised  classification tasks \cite{rasmus2015semi,zhao2015stacked}. However, \textit{DEPICT} is different from previous studies, since it is designed for the unsupervised clustering task, it also does not require max-pooling switches used in stacked what-where autoencoder (SWWAE) \cite{zhao2015stacked}, and lateral (skip) connections between encoder and decoder layers used in ladder network \cite{rasmus2015semi}. Algorithm \ref{alg:1} shows a brief description of \textit{DEPICT} algorithm.

\begin{algorithm}[!t]
	Initialize $\mathbf{Q}$ using a clustering algorithm\\
	\While{not converged  }{
		$\min\limits_{\boldsymbol{\psi}} \ \  -\frac{1}{N} \sum\limits_{ik} q_{ik} \log  \tilde{p}_{ik}  + \frac{1}{N}  \sum\limits_{il} \frac{1}{|\mathbf{z}_i^l|} \| \mathbf{z}_i^l - \hat{\mathbf{z}}_i^l \|_2^2$\\
		$ p_{ik}^{(t)} \propto exp( \boldsymbol{\theta}_k^T \mathbf{z}_i^{L} )$  \\
		$ q_{ik}^{(t)}   \propto   p_{ik}   / (\sum_{i'} p_{i'k})^{\frac{1}{2}}   $
	}
	\caption{\textit{DEPICT} Algorithm}
	\label{alg:1}
\end{algorithm}

\section{Experiments} \label{section:experiment}
In this section, we first evaluate \textit{DEPICT}\footnote{Our code is available in \url{https://github.com/herandy/DEPICT}} in comparison with state-of-the-art clustering methods on several benchmark image datasets. Then, the running speed of the best clustering models are compared. Moreover, we examine different learning approaches for training \textit{DEPICT}. Finally, we analyze the performance of \textit{DEPICT} model on semi-supervised classification tasks.

\noindent \textbf{Datasets}:
In order to show that \textit{DEPICT} works well with various kinds of datasets, we have chosen the following handwritten digit and face image datasets. Considering that clustering tasks are fully unsupervised, we concatenate the training and testing samples when applicable.
\textit{MNIST-full}: A dataset containing a total of 70,000 handwritten digits with 60,000 training and 10,000 testing samples, each being a 32 by 32 monochrome image \cite{lecun1998gradient}.
\textit{MNIST-test}: A dataset which only consists of the testing part of \textit{MNIST-full} data.
\textit{USPS}: It is a handwritten digits dataset from the \textit{USPS} postal service, containing 11,000 samples of 16 by 16 images.
\textit{CMU-PIE}: A dataset including 32 by 32 face images of 68 people with 4 different expressions \cite{sim2002cmu}.
\textit{Youtube-Face (YTF)}: Following \cite{yang2016joint}, we choose the first 41 subjects of YTF dataset. Faces inside images are first cropped and then resized to 55 by 55 sizes \cite{wolf2011face}.
\textit{FRGC}: Using the 20 random selected subjects in \cite{yang2016joint} from the original dataset, we collect 2,462 face images. Similarly, we first crop the face regions and resize them into 32 by 32 images.
Table \ref{table:datasets} provides a brief description of each dataset.

\begin{table}[h]
\centering
	\small
	\setlength{\tabcolsep}{7pt}
	\renewcommand{\arraystretch}{0.9}
	\begin{tabular}{l|c|c|c}
		\hline
		{Dataset}  & \# Samples 	& \# Classes 	& \# Dimensions \\ \hline \hline

		\textit{MNIST-full} 	& 70,000 		& 10		& 1$\times$28$\times$28 	  \\
		\textit{MNIST-test}	& 10,000 		& 10		& 1$\times$28$\times$28 	  \\
		\textit{USPS} 	& 11,000 		& 10		& 1$\times$16$\times$16 	  \\
		\textit{FRGC}	& 2,462 		& 20		& 3$\times$32$\times$32 	  \\
		\textit{YTF}	& 10,000 		& 41		& 3$\times$55$\times$55 	  \\
		\textit{CMU-PIE}	& 2,856 		& 68		& 1$\times$32$\times$32 	  \\

		\hline
	\end{tabular}

	\caption{Dataset Descriptions}
	\label{table:datasets}
\end{table}

\noindent \textbf{Clustering Metrics}:
We have used 2 of the most popular evaluation criteria widely used for clustering algorithms, accuracy (ACC) and normalized mutual information (NMI).
The best mapping between cluster assignments and true labels is computed using the Hungarian algorithm \cite{kuhn1955hungarian} to measure accuracy.
NMI calculates the normalized measure of similarity between two labels of the same data \cite{xu2003document}. Results of NMI do not change by permutations of clusters (classes), and they are normalized to have $[0,1]$ range, with $0$ meaning no correlation and $1$ exhibiting perfect correlation.

\begin{table*}[t]
	\small
	\setlength{\tabcolsep}{4pt}
	\renewcommand{\arraystretch}{0.95}
	\centering
	\begin{tabular}{l|cc|cc|cc|cc|cc|cc|c}
		\hline
		Dataset &
		\multicolumn{2}{c|}{\textit{MNIST-full}} &
		\multicolumn{2}{c|}{\textit{MNIST-test}} &
		\multicolumn{2}{c|}{\textit{USPS}} &
		\multicolumn{2}{c|}{\textit{FRGC}} &
		\multicolumn{2}{c|}{\textit{YTF}} &
		\multicolumn{2}{c|}{\textit{CMU-PIE}} &
		\multirow{2}{*}{\begin{tabular}{c}\# tuned \\ HPs\end{tabular}}\\ \cline{0-12}
		& NMI & ACC & NMI & ACC & NMI & ACC & NMI & ACC& NMI & ACC& NMI & ACC & \\
		\hline
		\textit{$K$-means} 	& 0.500$^{*}$ 	& 0.534$^{*}$ 	& 0.501$^{*}$ 	& 0.547$^{*}$ 	& 0.450$^{*}$ 	& 0.460$^{*}$ 	& 0.287$^{*}$ 	& 0.243$^{*}$ 	& 0.776$^{*}$ 	& 0.601$^{*}$ 	& 0.432$^{*}$ 	& 0.223$^{*}$ 	&	0	\\
		\textit{N-Cuts} 	& 0.411 		& 0.327 		& 0.753 		& 0.304 		& 0.675 		& 0.314 		& 0.285 		& 0.235 		& 0.742 		& 0.536 		& 0.411 		& 0.155 		&		0		\\
		\textit{SC-ST} 		& 0.416 		& 0.311 		& 0.756 		& 0.454 		& 0.726 		& 0.308 		& 0.431 		& 0.358			& 0.620			& 0.290 		& 0.581 		& 0.293 		&		0		\\
		\textit{SC-LS} 		& 0.706 		& 0.714		 	& 0.756 		& 0.740			& 0.681			& 0.659 		& 0.550 		& 0.407 		& 0.759 		& 0.544 		& 0.788 		& 0.549 		&		0		\\
		\textit{AC-GDL} 	& 0.017 		& 0.113 		& 0.844 		& 0.933			& 0.824			& 0.867 		& 0.351 		& 0.266 		& 0.622 		& 0.430			& 0.934 		& 0.842 		&		1		\\
		\textit{AC-PIC} 	& 0.017 		& 0.115 		& 0.853 		& 0.920			& 0.840			& 0.855 		& 0.415 		& 0.320 		& 0.697 		& 0.472 		& 0.902 		& 0.797 		&		0	\\
		\textit{SEC} 		& 0.779$^{*}$ 	& 0.804$^{*}$ 	& 	0.790$^{*}$	& 	0.815$^{*}$	& 0.511$^{*}$	& 0.544$^{*}$ 	&  	-			&  	-			&  	 -			&  	-			&	-			&	-			&	1	\\
		\textit{LDMGI} 		& 0.802$^{*}$ 	& 0.842$^{*}$ 	& 	0.811$^{*}$	& 	0.847$^{*}$	& 0.563$^{*}$	& 0.580$^{*}$	&  	-			&  	-			& 	 -			&	-			&	-			&	-			&	1	\\  \hline \hline
		\textit{NMF-D} 		& 0.152$^{*}$ 	& 0.175$^{*}$ 	& 0.241$^{*}$ 	& 0.250$^{*}$ 	& 0.287$^{*}$	& 0.382$^{*}$ 	& 0.259$^{*}$ 	& 0.274$^{*}$	& 0.562$^{*}$ 	& 0.536$^{*}$ 	& 0.920$^{*}$ 	& 0.810$^{*}$ 	& 0 \\
		\textit{TSC-D} 		& 0.651 		& 0.692			&  	-			& 	-			& 	-			&  	-			&   -			&  	-			&  	 -			& 	-			& 	-			&  - 			& 2	\\
		\textit{DEC} 		& 0.816$^{*}$ 	& 0.844$^{*}$ 	&  	0.827$^{*}$	& 	0.859$^{*}$	& 0.586$^{*}$ 	& 0.619$^{*}$	&  	0.505$^{*}$	&  	0.378$^{*}$	&  	 0.446$^{*}$	&  0.371$^{*}$	&	0.924$^{*}$	&	0.801$^{*}$ &	1	\\
		\textit{JULE-SF} 	& 0.906 		& 0.959 		& 0.876 		& 0.940			& 0.858			& 0.922 		& 0.566 		& 0.461 		 &\textbf{0.848} &\textbf{0.684} &0.984 			& 0.980 		&  3 \\
		\textit{JULE-RC} 	& 0.913 		& 0.964 		&\textbf{0.915} & 0.961 		& 0.913			& 0.950 		& 0.574 		& 0.461 		 &\textbf{0.848} &\textbf{0.684} & \textbf{1.00} & \textbf{1.00} & 3 \\ \hline \hline
		\textit{DEPICT}		&\textbf{0.917} &\textbf{0.965} &  \textbf{0.915} 		&\textbf{0.963} &\textbf{0.927} &\textbf{0.964} &\textbf{0.610} & \textbf{0.470}& 0.802 		& 0.621 		& 0.974 		& 0.883  		&	0 \\
		\hline
	\end{tabular}
	\vspace{-0.15cm}
	\caption{Clustering performance of different algorithms on image datasets based on accuracy (ACC) and normalized mutual information (NMI). The numbers of tuned hyper-parameters (\# tuned HPs) using the supervisory signals are also shown for each algorithm. The results of alternative models are reported from original papers, except the ones marked by (${*}$) on top, which are obtained by us running the released code. We put dash marks (-) for the results that are not practical to obtain. }
	\label{table:2}
	\vspace{-0.2cm}
\end{table*}

\subsection{Evaluation of Clustering Algorithm}

\noindent \textbf{Alternative Models}:
We compare our clustering model, \textit{DEPICT}, with several baseline and state-of-the-art clustering algorithms, including \textit{$K$-means}, normalized cuts (\textit{N-Cuts}) \cite{shi2000normalized}, self-tuning spectral clustering (\textit{SC-ST}) \cite{zelnik2004self}, large-scale spectral clustering (\textit{SC-LS}) \cite{chen2011large}, graph degree linkage-based agglomerative clustering (\textit{AC-GDL}) \cite{zhang2012graph}, agglomerative clustering via path integral (\textit{AC-PIC}) \cite{zhang2013agglomerative}, spectral embedded clustering (\textit{SEC}) \cite{nie2011spectral}, local discriminant models and global integration (\textit{LDMGI}) \cite{yang2010image}, \textit{NMF} with deep model (\textit{NMF-D}) \cite{trigeorgis2014deep}, task-specific clustering with deep model (\textit{TSC-D}) \cite{wang2016learning}, deep embedded clustering (\textit{DEC}) \cite{xie2016unsupervised}, and joint unsupervised learning (\textit{JULE}) \cite{yang2016joint}.

\noindent \textbf{Implementation Details}:
We use a common architecture for \textit{DEPICT} and avoid tuning any hyper-parameters using the labeled data in order to provide a practical algorithm for real-world clustering tasks. For all datasets, we consider two convolutional layers followed by a fully connected layer in encoder and decoder pathways. While for all convolutional layers, the feature map size is 50 and the kernel size is about $5 \times 5$, the dimension of the embedding subspace is set equal to the number of clusters in each dataset. We also pick the proper stride, padding and crop to have an output size of about $10 \times 10$ in the second convolutional layer. Inspired by \cite{radford2015unsupervised}, we consider leaky rectified (leaky RELU) non-linearity \cite{maas2013rectifier} as the activation function of convolutional and fully connected layers, except in the last layer of encoder and first layer of decoder, which have Tanh non-linearity functions. Consequently, we normalize the image intensities to be in the range of $[-1, 1]$.
Moreover, we set the learning rate and dropout to $10^{-4}$ and $0.1$ respectively, adopt adam as our optimization method with the default hyper-parameters $\beta_1 = 0.9$, $\beta_2 = 0.999$, $\epsilon=1e-08$ \cite{kingma2014adam}. The weights of convolutional and fully connected layers are all initialized by Xavier approach \cite{glorot2010understanding}.
Since the clustering assignments in the first iterations are random and not reliable for clustering loss, we first train \textit{DEPICT} without clustering loss function for a while, then initialize the clustering assignment $q_{ik}$ by clustering the embedding subspace features via simple algorithms like $K$-means or \textit{AC-PIC}.

\noindent \textbf{Quantitative Comparison}:
We run \textit{DEPICT} and other clustering methods on each dataset. We followed the implementation details for \textit{DEPICT} and report the average results from 5 runs. For the rest, we present the best reported results either from their original papers or from \cite{yang2016joint}. For unreported results on specific datasets, we run the released code with hyper-parameters mentioned in the original papers, these results are marked by (${*}$) on top.  But, when the code is not publicly available, or running the released code is not practical, we put dash marks (-) instead of the corresponding results.
Moreover, we mention the number of hyper-parameters that are tuned using supervisory signals (labeled data) for each algorithm. Note that this number only shows the quantity of hyper-parameters, which are set differently for various datasets for better performance.

Table \ref{table:2} reports the clustering metrics, normalized mutual information (NMI) and accuracy (ACC), of the algorithms on the aforementioned datasets. As shown, \textit{DEPICT} outperforms other algorithms on four datasets and achieves competitive results on the remaining two. It should be noted that we think hyper-parameter tuning using supervisory signals is not feasible in real-world clustering tasks, and hence \textit{DEPICT} is a significantly better clustering algorithm compared to the alternative models in practice. For example, \textit{DEC}, \textit{SEC}, and \textit{LDMGI} report their best results by tuning one hyper-parameter over nine different options, and \textit{JULE-SF} and \textit{JULE-RC} achieve their good performance by tweaking several hyper-parameters over various datasets. However, we do not tune any hyper-parameters for \textit{DEPICT} using the labeled data and only report the result with the same (default) hyper-parameters for all datasets.

\subsection{Running Time Comparison}

	In order to evaluate the efficiency of our clustering algorithm in dealing with large-scale and high dimensional data, we compare the running speed of \textit{DEPICT} with its competing algorithms, \textit{JULE-SF} and \textit{JULE-RC}.
	Moreover, the fast versions of \textit{\textit{JULE-SF}} and \textit{JULE-RC} are also evaluated. Note that JULE-SF(fast) and \textit{JULE-RC}(fast) both require tuning one extra hyper-parameter for each dataset to achieve results similar to the original \textit{JULE} algorithms in Table 2 \cite{yang2016joint}.
	 We run \textit{DEPICT} and the released code for \textit{JULE} algorithms\footnote{\url{https://github.com/jwyang/JULE-Torch}}  on a machine with one Titan X pascal GPU and a Xeon E5-2699 CPU.

	Figure \ref{fig:5} illustrates the running time for \textit{DEPICT} and \textit{JULE} algorithms on all datasets. Note that running times of \textit{JULE-SF} and \textit{JULE-RC} are shown linearly from 0 to 30,000 and logarithmically for larger values for the sake of readability. In total, \textit{JULE-RC}, \textit{JULE-SF}, \textit{JULE-RC}(fast), \textit{JULE-SF}(fast) and \textit{DEPICT} take 66.1, 35.5, 11.0, 6.6 and 4.7 hours respectively to run over all datasets.  While all algorithms have approximately similar running times on small datasets (\textit{FRGC} and \textit{CMU-PIE}), when dealing with the large-scale and high-dimensional datasets (\textit{MNIST-full} and \textit{YTF}), \textit{DEPICT} almost shows a linear increase in the running time, but the running times of original JULE algorithms dramatically grow with the size and number of input data.  This outcome again emphasizes the practicality of \textit{DEPICT} for real-world clustering tasks.
	\begin{figure}[t]
\centering
		\begin{center}
			\includegraphics[trim=0mm 6mm 19.5mm 12mm, clip,  width=0.90\linewidth]{figures/half_log_5-eps-converted-to.pdf}
		\end{center}
		\vspace{-0.5cm}
		\caption{Running time comparison of \textit{DEPICT} and \textit{JULE} clustering algorithms on image datasets.}
		\vspace{-0.4cm}
		\label{fig:5}
	\end{figure}

\begin{table*}[t]
	\small
	\setlength{\tabcolsep}{5pt}
	\renewcommand{\arraystretch}{0.95}
	\centering
	\begin{tabular}{ll|ll|ll|ll|ll|ll|ll}
		\hline
		Dataset &
		&\multicolumn{2}{c|}{\textit{MNIST-full}} &
		\multicolumn{2}{c|}{\textit{MNIST-test}} &
		\multicolumn{2}{c|}{\textit{USPS}} &
		\multicolumn{2}{c|}{\textit{FRGC}} &
		\multicolumn{2}{c|}{\textit{YTF}} &
		\multicolumn{2}{c}{\textit{CMU-PIE}}\\ \hline
		&& NMI & ACC & NMI & ACC & NMI & ACC & NMI & ACC& NMI & ACC & NMI & ACC\\
		\hline
		\multirow{3}{*}{\begin{minipage}{2.2cm} \centering  \textit{Deep-ConvAE} + \hspace{1.5cm} \textit{AC-PIC} \end{minipage} }
		&	SdA	&	0.255	&	0.348	&	0.313	&	0.345	&	0.223	&	0.290	&	0.120	&	0.230	&	0.414	&	0.302	&	0.354	&	0.266	 \\
		&	RdA	&	0.615	&	0.455	&	0.859	&	0.900	&	0.886	&	0.866	&	0.443	&	0.363	&	0.597	&	0.425	&	0.912	&	0.817	 \\
		&	MdA	&	0.729	&	0.506	&	0.876	&	0.942	&	0.906	&	0.878	&	0.583	&	0.427	&	0.640	&	0.448	&	0.931	&	0.883	 \\
		\hline
		\hline
		\multirow{3}{*}{ \centering \textit{DEPICT}}
		&	SdA	&	0.365	&	0.427	&	0.353	&	0.390	&	0.328	&	0.412	&	0.211	&	0.300	&	0.414	&	0.302	&	0.354	&	0.266	 \\
		&	RdA	&	0.808	&	0.677	&	0.899	&	0.950	&	0.901	&	0.923	&	0.551	&	0.444	&	0.652	&	0.450	&	0.951	&	0.926	 \\
		& MdA  & \textbf{0.917} & \textbf{0.965} &  \textbf{0.915} & \textbf{0.963} & \textbf{0.927} & \textbf{0.964}  & \textbf{0.610} & \textbf{0.470} & \textbf{0.802} & \textbf{0.621} & \textbf{0.974} & \textbf{0.883}  \\
		\hline
	\end{tabular}
	\vspace{-0.25cm}
	\caption{Clustering performance of different learning approaches, including SdA, RdA and MdA, for training \textit{DEPICT} and \textit{Deep-ConvAE}+\textit{AC-PIC} models.}
	\label{table:3}
	\vspace{-0.35cm}
\end{table*}

\subsection{Evaluation of Learning Approach}
In order to evaluate our joint learning approach, we compare  several strategies for training \textit{DEPICT}. For training a multi-layer convolutional autoencoder, we analyze the following three approaches : 1) Standard stacked denoising autoencoder (SdA), in which the model is first pretrained using the reconstruction loss function in a layer-wise manner, and the encoder pathway is then fine-tuned using the clustering objective function \cite{vincent2010stacked}. 2) Another approach (RdA)  is suggested in \cite{xie2016unsupervised} to improve the SdA learning approach, in which all of the autoencoder layers are retrained after the pretraining step, only using the reconstruction of input layer while data is not corrupted by noise. The fine-tuning step is also done after the retraining step. 3) Our learning approach (MdA), in which the whole model is trained simultaneously using the joint reconstruction loss functions from all layers along with the clustering objective function.

Furthermore, we also examine the effect of clustering loss (through error back-prop) in constructing the embedding subspace. To do so, we train a similar multi-layer convolutional autoencoder (\textit{Deep-ConvAE}) only using the reconstruction loss function to generate the embedding subspace. Then, we run the best shallow clustering algorithm (\textit{AC-PIC}) on the embedded data. Hence, this model (\textit{Deep-ConvAE}+\textit{AC-PIC}) differs from \textit{DEPICT} in the sense that its embedding subspace is only constructed using the reconstruction loss and does not involve the clustering loss.

Table \ref{table:3} indicates the results of \textit{DEPICT} and \textit{Deep-ConvAE}+\textit{\textit{AC-PIC}} when using the different learning approaches. As expected, \textit{DEPICT} trained by our joint learning approach (MdA) consistently outperforms the other alternatives on all datasets. Interestingly, MdA learning approach shows promising results for \textit{Deep-ConvAE}+\textit{AC-PIC} model, where only reconstruction losses are used to train the embedding subspace. Thus, our learning approach is an efficient strategy for training autoencoder models due to its superior results and fast end-to-end training.

\subsection{Semi-Supervised Classification Performance}

	Representation learning in an unsupervised manner or using a small number of labeled data has recently attracted great attention. Due to the potential of our model in learning a discriminative embedding subspace, we evaluate \textit{DEPICT} in a semi-supervised classification task. Following the semi-supervised experiment settings \cite{rasmus2015semi,zhao2015stacked}, we train our model using a small random subset of  \textit{MNIST-training} dataset as labeled data and the remaining as unlabeled data. The classification error of \textit{DEPICT} is then computed using the \textit{MNIST-test} dataset, which is not seen during training.  Compared to our unsupervised learning approach, we only utilize the clusters corresponding to each labeled data in training process. In particular, only for labeled data, the cluster labels (assignments) are set using the best map technique from the original classification labels once, and then they will be fixed during the training step.

	Table \ref{table:5} shows the error results for several semi-supervised classification models using different numbers of labeled data. Surprisingly, \textit{DEPICT} achieves comparable results with the state-of-the-art, despite the fact that the semi-supervised classification models use 10,000 validation data to tune their hyper-parameters, \textit{DEPICT} only employs the labeled training data (e.g. 100) and does not tune any hyper-parameters. Although \textit{DEPICT} is not mainly designed for classification tasks, it outperforms several models including \textit{SWWAE} \cite{zhao2015stacked}, \textit{M1+M2} \cite{kingma2014semi}, and \textit{AtlasRBF} \cite{pitelis2014semi}, and has comparable results with the complicated \textit{Ladder} network \cite{rasmus2015semi}. These results further confirm the discriminative quality of the embedding features of \textit{DEPICT}.

\begin{table}[t]
	\small
	\setlength{\tabcolsep}{6pt}
	\renewcommand{\arraystretch}{0.90}
	\centering
	\begin{tabular}{l|ccc}
		\hline
		Model							& 100			& 1000 			& 3000 	\\\hline
		\textit{T-SVM} \cite{vapnik1998statistical}	& 16.81 		& 5.38			& 3.45	\\
		\textit{CAE} \cite{rifai2011contractive}		& 13.47			& 4.77			& 3.22	\\
		\textit{MTC} \cite{rifai2011manifold}		& 12.03 		& 3.64  		& 2.57	\\
		\textit{PL-DAE} \cite{lee2013pseudo}			& 10.49			& 3.46			& 2.69\\
		\textit{AtlasRBF} \cite{pitelis2014semi}		& 8.10			& 3.68			& - \\
		\textit{M1+M2} \cite{kingma2014semi}			& 3.33$\pm$0.14	& 2.40$\pm$0.05	& 2.18$\pm$0.04	\\
		\textit{SWWAE} \cite{zhao2015stacked}		& 8.71$\pm$0.34	& 2.83$\pm$0.10	& 2.10$\pm$0.22	\\
		\textit{Ladder} \cite{rasmus2015semi}		& \textbf{1.06$\pm$0.37}	& \textbf{0.84$\pm$0.08}	& 	-	\\ \hline
		\textit{DEPICT} 								& 2.65$\pm$0.35	& 2.10$\pm$0.11	& \textbf{1.91$\pm$0.06}	\\ \hline
	\end{tabular}
	\vspace{-0.15cm}
	\caption{Comparison of \textit{DEPICT} and several semi-supervised classification models in \textit{MNIST} dataset with different numbers of labeled data. }
	\label{table:5}
	\vspace{-0.4cm}
\end{table}

\section{Conclusion}\label{section:conclusion}
In this paper, we proposed a new deep clustering model, \textit{DEPICT}, consisting of a soft-max layer stacked on top of a multi-layer convolutional autoencoder. We employed a regularized relative entropy loss function for clustering, which leads to balanced cluster assignments. Adopting our autoencoder reconstruction loss function enhanced the embedding learning. Furthermore, a joint learning framework was introduced to train all network layers simultaneously and avoid layer-wise pretraining. Experimental results showed that \textit{DEPICT} is a good candidate for real-world clustering tasks, since it achieved superior or competitive results compared to alternative methods while having faster running speed and not needing hyper-parameter tuning. Efficiency of our joint learning approach was also confirmed in clustering and semi-supervised classification tasks.

{\small
	\bibliographystyle{ieee}


text:
\section{Introduction} \label{sec:intr}
There has been vast research on the non-vanishing of central $L$-values for families of automorphic forms, since the pioneering work of Duke \cite{duke1995} and Iwaniec--Sarnak \cite{iwaniec-sarnak1999dirichlet,iwaniec-sarnak2000siegel}. To get \textit{positive-proportional} non-vanishing results in families, one typically turns to the method of moments and the mollification method \`{a} la Selberg (see, for example, \cite{KowalskiMichel1999, VanderKam1999, KowalskiMichel2000, KMV2000, KMV2000b, Soundararajan2000, KMV2002, Blomer2008, Khan2010, Luo2015}, and others). In the current work we follow this approach and go beyond families of GL(2) forms (and symmetric-square lifts of GL(2) forms), and study the central $L$-values of Maass forms on GL(3) and prove a non-vanishing result of such values (Theorem \ref{main result}), which is a positive-proportional result in the sense of Remark \ref{remark_main_result}.

To state our result, we introduce a few notations and refer the reader to \S\,\ref{subsec:cusp} for certain details. Pick an orthogonal basis $\{\phi_j\}$ of Hecke--Maass forms for $\Gamma=SL(3,\mathbb{Z})$. Each $\phi_j$ has spectral parameter $\nu_{j}=\big(\nu_{j,1},\nu_{j,2},\nu_{j,3}\big)$, the Langlands parameter
$\mu_{j}=\big(\mu_{j,1},\mu_{j,2},\mu_{j,3}\big)$, and the Hecke eigenvalues $A_j(n,1)$. The main objects under investigation are the $L$-functions
\[
  L(s,\phi_j) = \sum_{n=1}^{\infty} \frac{A_{j}(1,n)}{n^{s}}

    \quad \textrm{for } \Re(s)>1.
\]
A simple observation is that there is no trivial reason for $L(\frac{1}{2},\phi_j)$ to vanish, since every $\phi_j$ is necessarily even and the sign of the functional equation of $L(s,\phi_j)$ is positive. In fact, one expects many of $L(\frac{1}{2},\phi_j)$ to be nonzero.
As in Blomer--Buttcane \cite{blomer2015subconvexity}, we consider the generic case in short interval.
Let $\mu_0=(\mu_{0,1},\mu_{0,2},\mu_{0,3})$ and $\nu_0=(\nu_{0,1},\nu_{0,2},\nu_{0,3})$, and satisfy
the corresponding relations \eqref{eqn:m2n} and \eqref{eqn:n2m}.
We also assume
\[
  |\mu_{0,i}|\asymp|\nu_{0,i}|\asymp T:=\|\mu_0\|, \quad
  1\leq i \leq 3.
\]
Let $M=T^\theta$ for any fixed $0<\theta<1$. Define a test function $h_{T,M}(\mu)$ (depending on $\mu_0$) for $\mu=(\mu_1,\mu_2,\mu_3)\in\mathbb{C}^3$ by
$$
     h_{T,M}(\mu) := P(\mu)^2 \bigg(\sum_{w\in\cW}\psi\bigg(\frac{w(\mu)-\mu_0}{M}\bigg)\bigg)^2,
$$
where
$$
    \psi(\mu) =\exp\left(\mu_{1}^2+\mu_{2}^2+\mu_{3}^2\right)
$$
and
$$
    P(\mu)=\prod_{0\leq n\leq A} \prod_{k=1}^{3} \frac{\left(\nu_{k}-\frac{1}{3}(1+2n)\right)\left(\nu_{k}+\frac{1}{3}(1+2n)\right)}{|\nu_{0,k}|^2}
$$
for some fixed large $A>0$. Here
\[
  \mathcal{W} := \left\{ I,\; w_2=\left(\begin{smallmatrix} 1 & & \\  & & 1 \\  &1&   \end{smallmatrix}\right),\;
  w_3=\left(\begin{smallmatrix}  &1& \\  1&& \\  &&1   \end{smallmatrix}\right),\;
  w_4=\left(\begin{smallmatrix}  &1& \\  &&1 \\  1&&   \end{smallmatrix}\right),\;
  w_5=\left(\begin{smallmatrix}  &&1 \\  1&& \\  &1&   \end{smallmatrix}\right),\;
  w_6=\left(\begin{smallmatrix}  &&1 \\  &1& \\  1&&   \end{smallmatrix}\right)
  \right\}
\]
is the Weyl group of $SL(3,\bbR)$. The function $h(\mu)$ has the localizing effect at a ball of radius $M$ about $w(\mu_0)$ for each $w\in\cW$, and other nice properties stated in \S\,\ref{subsec:cusp}. Then with the normalizing factor
$$
    \cN_j:=\|\phi_j\|^2\prod_{k=1}^{3}\cos\left(\frac{3}{2}\pi\nu_{j,k}\right),
$$
our main result is as follows.
\begin{theorem}\label{main result}
  We have
  \[
    \sum_{\substack{j\geq1 \\ L(\thf,\phi_j)\neq0}} \frac{h_{T,M}(\mu_j)}{\cN_j} \gg T^3M^2.
  \]
\end{theorem}

\begin{remark}\label{remark_main_result}
  By a stronger form of the GL(3) spectral large sieve inequality obtained by Young (\cite[Theorem 1.1]{young2016largesieve}), one can get the following weighted Weyl law:
  \[
    \sum_{j\geq1} \frac{h_{T,M}(\mu_j)}{\cN_j} \ll T^3M^2.
  \]
  In fact, we can replace the above ``$\ll$" by ``$\asymp$".
  Thus in this sense Theorem \ref{main result} gives a positive-proportional non-vanishing result in short interval.
\end{remark}

Next we outline the structure of the paper and give the proof of Theorem \ref{main result}. In \S\,\ref{subsec:cusp} we briefly review facts of Maass forms and their $L$-functions, as well as the main analytic tool, the GL(3) Kuznetsov trace formula (Lemma \ref{lemma: KTF}). Define the mollifier $M_j$ for $L(\hf,\phi_j)$ by
\[
  M_j := \sum_{\ell\leq L}\frac{A_j(1,\ell)}{\ell^{1/2}} x_{\ell},
\]
where
\begin{align}
  x_{\ell} := \mu(\ell)\f{1}{2\pi i}\int_{(3)}\f{(L/\ell)^s}{s^2}\frac{\dd s}{\log L}
  =  \begin{cases}
        \mu(\ell)\f{\log(L/\ell)}{\log L},\ \ &\ell\leq L,
        \\\noalign{\vskip 1mm}
        0,\ \ \ &\ell>L,
     \end{cases}
\end{align}
and $L=T^{\delta}$ for some small $\delta>0$. Then we study the mollified moments of the central $L$-values and prove the following two propositions
in \S\ref{sec:2ndmoment} and \S\ref{sec:1stmoment}, respectively.

\begin{proposition}\label{second moment}
  We have
  \[
    \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} |L(\thf,\phi_j) M_j|^2 \ll T^3M^2,
  \]
  provided $0<\delta<11/78$.
\end{proposition}

\begin{proposition}\label{first moment}
  We have
  \[
    \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} L(\thf,\phi_j) M_j \asymp T^3M^2,
  \]
  provided $0<\delta<11/78$.
\end{proposition}
\begin{remark}
  The above results can be improved. The restriction of $\delta$ comes from the contribution of Eisenstein series, which can be refined if we use the subconvexity bounds for GL(1) and GL(2) $L$-functions or the average Lindel\"{o}f bound of the related families of $L$-functions.
\end{remark}

Through out the paper, $\varepsilon$ is an arbitrarily small positive number
and $B$ is a sufficiently large positive number which may not be the same at each occurrence.

\section{Preliminaries}
In this section we review essential facts and tools, required for later development.
\subsection{Hecke--Maass cusp forms and their $L$-functions}\label{subsec:cusp}
Let $G=GL(3,\mathbb{R})$ with maximal compact subgroup $K=O(3,\mathbb{R})$ and center $Z\cong\mathbb{R}^\times$. Let $\mathbb{H}_3=G/(K\cdot Z)$ be the generalized upper half-plane.
For $0\leq c\leq \infty$, let
\begin{equation*}
  \Lambda_{c} = \Big\{ \mu=(\mu_1,\mu_2,\mu_3)\in\mathbb{C}^3 \; \Big| \; \mu_1+\mu_2+\mu_3=0,\ |\Re\mu_k|\leq c,\ k=1,2,3 \Big\},
\end{equation*}
and
\begin{equation*}
  \Lambda_{c}' = \Big\{ \mu\in\Lambda_c \; \Big| \; \{-\mu_1,-\mu_2,-\mu_3\}=\{\overline{\mu_1},\overline{\mu_2},\overline{\mu_3}\} \Big\}.
\end{equation*}
Consider a Hecke--Maass form $\phi$ in $L^2(SL(3,\mathbb{Z})\backslash\mathbb{H})$. Here $\mu$ will be the Langlands parameters.
Define the spectral parameters
\begin{equation}\label{eqn:n2m}
  \nu_1=\frac{1}{3}(\mu_1-\mu_2),\quad
  \nu_2=\frac{1}{3}(\mu_2-\mu_3),\quad
  \nu_3=\frac{1}{3}(\mu_3-\mu_1).
\end{equation}
We have
\begin{equation}\label{eqn:m2n}
  \mu_1=2\nu_1+\nu_2, \quad
  \mu_2=\nu_2-\nu_1, \quad
  \mu_3=-\nu_1-2\nu_2.
\end{equation}
We will simultaneously use $\mu$ and $\nu$.
By unitarity and the standard Jacquet--Shalika bounds, the Langlands parameter of an
arbitrary irreducible representation $\pi\subseteq L^2(SL(3,\mathbb{Z})\backslash\mathbb{H}_3)$
is contained in $\Lambda_{1/2}'\subseteq \Lambda_{1/2}$,
and the non-exceptional parameters are in $\Lambda_{0}'=\Lambda_{0}$.

Let $\phi$ be a Hecke--Maass cusp form for $\Gamma=SL(3,\mathbb{Z})$ with Fourier coefficients
$A_{\phi}(m,n)\in\mathbb{C}$ for $(m,n)\in \mathbb{N}^2$.
The standard L-function of $\phi$ is given by
\[
  L(s,\phi):=\sum_{n=1}^{\infty} \frac{A_\phi(1,n)}{n^s} \quad \textrm{for }\Re(s)>1.
\]
For the dual form $\widetilde{\phi}$ the coefficients of $L(s,\widetilde{\phi})$ are $A_\phi(1,n)=\overline{A_\phi(n,1)}$.
The functional equation of $L(s,\phi)$ is
\begin{equation}\label{eqn:FE}
  L(s,\phi)\prod_{j=1}^{3}\Gamma_{\mathbb{R}}(s-\mu_j)
  = L(1-s,\widetilde{\phi})\prod_{j=1}^{3}\Gamma_{\mathbb{R}}(1-s+\mu_j) ,
\end{equation}
where $\Gamma_{\mathbb{R}}(s) = \pi^{-s/2}\Gamma(s/2)$.

\subsection{The minimal Eisenstein series and its Fourier coefficients}

Let
$$
U_3=\left\{\begin{pmatrix} 1&*&* \\ 0&1&* \\ 0&0&1 \end{pmatrix}\right\} \cap \Gamma.
$$
For $z\in \mathbb{H}_3$ and $\Re(\nu_1),\Re(\nu_2)$ sufficiently large,
we define the minimal Eisenstein series
\[
  E(z;\mu_1,\mu_2) := \sum_{\gamma\in U_3\backslash\Gamma} I_{\nu_1,\nu_2}(\gamma z),
\]
where
\[
  I_{\nu_1,\nu_2}(z) := y_1^{1+\nu_1+2\nu_2}y_2^{1+2\nu_1+\nu_2},
\]
and
\begin{equation}\label{z}
  z=\left\{\begin{pmatrix} y_1y_2&y_1x_2&x_3 \\ 0&y_1&x_1 \\ 0&0&1 \end{pmatrix}\right\}
  =\left\{\begin{pmatrix} 1&x_2&x_3 \\ 0&1&x_1 \\ 0&0&1 \end{pmatrix}\right\}\left\{\begin{pmatrix} y_1y_2&0&0 \\ 0&y_1&0 \\ 0&0&1 \end{pmatrix}\right\},
\end{equation}
with $y_1,y_2>0$ and $x_1,x_2,x_3\in\mathbb{R}$.
It has meromorphic continuation in $\nu_1$ and $\nu_2$.
The Fourier coefficients $A_{\mu}(m_1,m_2)$ is defined by
(see Goldfeld \cite[Theorems 10.8.6]{goldfeld2006automorphic})
\[
  A_{\mu}(m,1) = \sum_{d_1d_2d_3=m} d_1^{\mu_1}d_2^{\mu_2}d_3^{\mu_3},
\]
and the symmetry and Hecke relation (see Goldfeld \cite[Theorems 6.4.11]{goldfeld2006automorphic})
\[
  \begin{split}
     A_{\mu_1,\mu_2}(m,1) & = \overline{A_{\mu_1,\mu_2}(1,m)}, \\
     A_{\mu_1,\mu_2}(m,n) & = \sum_{d|(m,n)} \mu(d) A_{\mu_1,\mu_2}(m/d,1)A_{\mu_1,\mu_2}(1,n/d).
  \end{split}
\]
Hence we have
\begin{equation}\label{eqn:Amin<<}
  |A_{\mu}(m_1,m_2)| \ll_{\varepsilon} (m_1m_2)^{\varepsilon},
  \quad \textrm{if}\
  \mu\in (i\mathbb{R})^2.
\end{equation}
In order to state the Kuznetsov trace formula in the \S 2.5, we introduce
$$
\mathcal{N}_\mu=\frac{1}{16}\prod_{k=1}^3|\zeta(1+3\nu_{k})|^2
$$
corresponding to the minimal Eisenstein series $E(z;\mu_1,\mu_2)$,
where $\mu=(\mu_1,\mu_2,\mu_3)$.
Recall that (see \cite{titchmarsh1986theory}) we have
$$
\frac{1}{\zeta(1+it)}\ll \log (1+|t|),
$$
which implies that
\begin{align}\label{N min}
  \frac{1}{\mathcal{N}_\mu}\ll \prod_{k=1}^3\log^2(1+|\nu_k|).
\end{align}

\subsection{The maximal Eisenstein series and its Fourier coefficients}

Let
$$
P_{2,1}=\left\{\begin{pmatrix} *&*&* \\ *&*&* \\ 0&0&* \end{pmatrix}\right\} \cap \Gamma.
$$
Let $\mu\in\mathbb{C}$ have sufficiently large real part,
and let $g: SL(2,\mathbb{Z})\backslash\mathbb{H}_2\rightarrow \mathbb{C}$
be a Hecke--Maass cusp form with $\|g\|=1$,
Langlands parameter $\mu_g\in i\mathbb{R}$
and Hecke eigenvalue $\lambda_g(m)$.
The maximal Eisenstein series twisted by a Maass form $g$ is defined by
\[
  E(z;\mu;g) := \sum_{\gamma \in P_{2,1}\backslash\Gamma}\det(\gamma z)^{1/2+\mu} g(\mathfrak{m}_{P_{2,1}}(\gamma z)),
\]
where $z$ is defined as in \eqref{z}, and
\[
  \mathfrak{m}_{P_{2,1}}: \mathbb{H}_3 \rightarrow \mathbb{H}_2,\quad
  \left(\begin{matrix}
          y_1y_2 & y_1x_2 & x_3 \\
           & y_1 & x_1 \\
           &   & 1
        \end{matrix}\right) \mapsto
  \left(\begin{matrix}
          y_2 & x_2  \\
           &   1
        \end{matrix}\right)
\]
is the restriction to the upper left corner.
It has a meromorphic continuation in $\mu$.
The Fourier coefficients are determined by
\[
  B_{\mu,g}(m,1) = \sum_{d_1d_2=m}\lambda_g(d_1)d_1^{\mu}d_2^{-2\mu},
\]
and the symmetry and Hecke relation as above
(see Goldfeld \cite[Proposition 10.9.3 and Theorem 6.4.11]{goldfeld2006automorphic}).
Recall that we have the following Kim--Sarnak bound for $GL(2)$ Fourier coefficients
(see Kim \cite[Appendix 2]{kim2003functoriality})
\[
  |\lambda_g(n)|\ll_{\varepsilon} n^{7/64+\varepsilon}.
\]
Hence we have
\begin{equation}\label{eqn:Amax<<}
  |B_{\mu,g}(m_1,m_2)| \ll_\varepsilon (m_1m_2)^{7/64+\varepsilon},
  \quad \textrm{if}\
  \mu\in i\mathbb{R}.
\end{equation}
We also introduce
$$
\mathcal{N}_{\mu,g}=8L(1,\textup{Ad}^2g)|L(1+3\mu,g)|^2,
$$
where $L(s,\textup{Ad}^2g)$ is the adjoint square $L$-function of $g$,
and $L(s,g)$ is the $L$-function of $g$.
We have the lower bounds
$$
L(1,\textup{Ad}^2g)\gg (1+|\mu_g|)^{-\varepsilon},  \quad
L(1+it,g)\gg (1+|t|+|\mu_g|)^{-\varepsilon}.
$$
These lower bounds follow from \cite{hoffstein-lockhart, hoffstein-ramakrishnan, jacquet-shalika}, and \cite{gelbart-lapid-sarnak}.
Therefore, for $\mu\in i\mathbb{R}$, it follows that
\begin{align}\label{N max}
  \frac{1}{\mathcal{N}_{\mu,g}}\ll (1+|\mu|+|\mu_g|)^{\varepsilon}.
\end{align}

\subsection{The Kloosterman sums}

For $n_1$, $n_2$, $m_1$, $m_2$, $D_1$, $D_2\in\mathbb{N}$, we need the relevant Kloosterman sums
\begin{align*}
  \tilde{S}(n_1,n_2,m_1;D_1,D_2) :=\sum_{\substack{C_1\,(\mod D_1),\, C_2\,(\mod D_2)\\ (C_1,D_1)=(C_2,D_2/D_1)=1}}
  e\left(n_2\f{\bar{C_1}C_2}{D_1}+m_1\f{\bar{C_2}}{D_2/D_1}+n_1\f{C_1}{D_1}\right)
\end{align*}
for $D_1|D_2$, and
\begin{align*}
  &S(n_1,m_2,m_1,n_2;D_1,D_2)
  \\
  &:=\sum_{\substack{B_1,C_1(\mod D_1)\\ B_2,C_2(\mod D_2)\\ D_1C_2+B_1B_2+D_2C_1\equiv0\,(\mod D_1D_2)\\ (B_j,C_j,D_j)=1}}e\left(\f{n_1B_1+m_1(Y_1D_2-Z_1B_2)}{D_1}+\f{m_2B_2+n_2(Y_2D_1-Z_2B_1)}{D_2}\right),
\end{align*}
where $Y_jB_j+Z_jC_j\equiv1\,(\mod D_j)$ for $j=1,2$.

\subsection{The Kuznetsov trace formula}

We first introduce some notation.
Define the spectral measure on the hyperplane $\mu_1+\mu_2+\mu_3=0$ by
\[
  \dd_{\spec}\mu=\spec(\mu)\dd\mu,
\]
where
$$
  \spec(\mu):= \prod_{k=1}^{3}\left(3\nu_k\tan\left(\frac{3\pi}{2}\nu_k\right)\right)\quad\mbox{and}\quad \dd\mu=\dd\mu_1\dd\mu_2=\dd\mu_2\dd\mu_3=\dd\mu_1\dd\mu_3.
$$
Following \cite[Theorems 2 \& 3]{buttcane2016spectral},
we define the following integral kernels in terms of
Mellin--Barnes representations.
For $s \in \bbC$, $\mu \in \Lambda_{\infty}$ define the meromorphic function
$$
  \tilde{G}^{\pm}(s, \mu) := \frac{  \pi^{-3s}}{12288 \pi^{7/2}}\Biggl(\prod_{k=1}^3\frac{\Gamma(\frac{1}{2}(s-\mu_k))}{\Gamma(\frac{1}{2}(1-s+\mu_k))} \pm i   \prod_{k=1}^3\frac{\Gamma(\frac{1}{2}(1+s-\mu_k))}{\Gamma(\frac{1}{2}(2-s+\mu_k))} \Biggr),
$$
and for $s = (s_1, s_2) \in \bbC^2$, $\mu \in \Lambda_{\infty}$ define the meromorphic function
$$
  G(s, \mu) :=  \frac{1}{\Gamma(s_1 + s_2)} \prod_{k=1}^3 \Gamma(s_1 - \mu_k) \Gamma(s_2 + \mu_k).
$$
The latter is essentially the double Mellin transform
of the GL(3) Whittaker function.
We also define the following trigonometric functions

\begin{displaymath}
  \begin{split}
    & S^{++}(s, \mu) := \frac{1}{24 \pi^2} \prod_{k=1}^3 \cos\left(\frac{3}{2} \pi \nu_k\right),\\
    &  S^{+-}(s, \mu) :=  -\frac{1}{32 \pi^2} \frac{\cos(\frac{3}{2} \pi \nu_2)\sin(\pi(s_1 - \mu_1))\sin(\pi(s_2 + \mu_2))\sin(\pi(s_2 + \mu_3))}{\sin(\frac{3}{2} \pi \nu_1)\sin(\frac{3}{2} \pi \nu_3) \sin(\pi(s_1+s_2))}, \\
    & S^{-+}(s, \mu) :=-\frac{1}{32 \pi^2}  \frac{\cos(\frac{3}{2} \pi \nu_1)\sin(\pi(s_1 - \mu_1))\sin(\pi(s_1 - \mu_2))\sin(\pi(s_2 + \mu_3))}{\sin(\frac{3}{2} \pi \nu_2)\sin(\frac{3}{2} \pi \nu_3)\sin(\pi(s_1+s_2))}, \\
    & S^{--}(s, \mu) := \frac{1}{32 \pi^2}  \frac{\cos(\frac{3}{2} \pi \nu_3) \sin(\pi(s_1 - \mu_2))\sin(\pi(s_2 + \mu_2))}{\sin(\frac{3}{2} \pi \nu_2)\sin(\frac{3}{2} \pi \nu_1)}.
  \end{split}
\end{displaymath}
For  $y \in \Bbb{R} \setminus \{0\}$ with
$\sgn(y) = \epsilon$,
let
$$
  K_{w_4}(y;\mu) :=  \int_{-i\infty}^{i\infty} | y|^{-s} \tilde{G}^{\epsilon}(s, \mu) \frac{\dd s}{2\pi i}.
$$
For $y = (y_1, y_2) \in (\bbR\setminus \{0\})^2$ with
$\sgn(y_1) = \epsilon_1$, $\sgn(y_2) = \epsilon_2$, let
\begin{displaymath}
  \begin{split}
     K^{\epsilon_1, \epsilon_2}_{w_6}(y; \mu)
     :=  & \int_{-i\infty}^{i\infty} \int_{-i\infty}^{i\infty}  |4\pi^2 y_1|^{-s_1} |4\pi^2 y_2|^{-s_2}  G(s, \mu) S^{\epsilon_1, \epsilon_2}(s, \mu)\frac{\dd s_1 \dd s_2}{(2\pi i)^2}.
  \end{split}
\end{displaymath}

We can now state the Kuznetsov trace formula in the version of
Buttcane \cite[Theorems 2, 3, 4]{buttcane2016spectral}.

\begin{lemma}\label{lemma: KTF}
  Let $n_1$, $n_2$, $m_1$, $m_2 \in \bbN$ and
  let $h$ be a function  that is holomorphic on
  $\Lambda_{1/2 + \delta}$ for some $\delta > 0$,
  symmetric under the Weyl group $\cW$, of rapid decay when
  $|\Im \mu_j| \rightarrow \infty$,  and satisfies
  \begin{equation}\label{eqn: zeros}
    h(3\nu_j \pm 1)  = 0, \quad j = 1, 2, 3.
  \end{equation}
  Then we have
  \begin{displaymath}
    \begin{split}
      &  \cC + \cE_{\min} + \cE_{\max}
       = \Delta + \Sigma_4 + \Sigma_5 + \Sigma_6,
    \end{split}
  \end{displaymath}
  where
  \begin{equation*}
    \begin{split}
       \cC & := \sum_{j} \frac{h(\mu_{j})}{\cN_j} \overline{A_{j}(m_1, m_2)} A_{j}(n_1, n_2), \\
        \cE_{\min} & := \frac{1}{24(2\pi i)^2} \iint_{\Re(\mu)=0} \frac{h(\mu)}{\cN_{\mu}}
                        \overline{A_{\mu}(m_1, m_2)} A_{\mu}(n_1, n_2) \dd\mu_1 \dd\mu_2,\\
       \cE_{\max} & :=  \sum_{g} \frac{1}{2\pi i} \int_{\Re(\mu)=0}
            \frac{h(\mu+\mu_g,\mu-\mu_g,-2\mu)}{\cN_{\mu,g}}
             \overline{B_{\mu,g}(m_1, m_2)} B_{\mu,g}(n_1, n_2) \dd\mu,\\
    \end{split}
  \end{equation*}
  and
  \begin{equation*
    \begin{split}
      \Delta& := \delta_{n_1, m_1} \delta_{n_2, m_2}  \frac{1}{192\pi^5} \int_{\Re \mu = 0} h(\mu) \dd_{\spec}\mu,\\
       \Sigma_{4}& := \sum_{\epsilon  = \pm 1} \sum_{\substack{D_2 \mid D_1\\  m_2 D_1= n_1 D_2^2}}\frac{ \tilde{S}(-\epsilon n_2, m_2, m_1; D_2, D_1)}{D_1D_2} \Phi_{w_4}\!\left(  \frac{\epsilon m_1m_2n_2}{D_1 D_2} \right),  \\
       \Sigma_{5} &:= \sum_{\epsilon  = \pm 1} \sum_{\substack{     D_1 \mid D_2\\ m_1 D_2 = n_2 D_1^2}} \frac{ \tilde{S}(\epsilon n_1, m_1, m_2; D_1, D_2) }{D_1D_2}\Phi_{w_5}\!\left( \frac{\epsilon n_1m_1m_2}{D_1 D_2}\right),\\
       \Sigma_6 &:= \sum_{\epsilon_1, \epsilon_2 = \pm 1} \sum_{D_1,  D_2  } \frac{S(\epsilon_2 n_2, \epsilon_1 n_1, m_1, m_2; D_1, D_2)}{D_1D_2} \Phi_{w_6}\!  \left( - \frac{\epsilon_2 m_1n_2D_2}{D_1^2}, - \frac{\epsilon_1 m_2n_1D_1}{ D_2^2}\right),
    \end{split}
  \end{equation*}
  with
\begin{equation}\label{defPhi}
    \begin{split}
      & \Phi_{w_4}(y) :=  \int_{\Re \mu = 0} h(\mu) K_{w_4}(y; \mu ) \dd_{\spec}\mu,\\
      & \Phi_{w_5}(y) := \int_{\Re \mu = 0} h(\mu) K_{w_4}(-y; -\mu ) \dd_{\spec}\mu,\\
      & \Phi_{w_6}(y_1, y_2) := \int_{\Re \mu = 0} h(\mu) K^{\sgn(y_1), \sgn(y_2)}_{w_6}((y_1, y_2) ;  \mu ) \dd_{\spec}\mu.
    \end{split}
\end{equation}
\end{lemma}

In the first moment, we will use $h(\mu)=h_{T,M}(\mu)$ to be the test function which is same as Blomer--Buttcane \cite{blomer2015subconvexity}, and in the second moment, we will use
$h_2(\mu)=h_{T,M}(\mu)W_{\mu,N}(m_1m_2)$ (see \eqref{test h2}) to be the test function.
The function $h(\mu)$ localizes at a ball of radius $M$ about $w(\mu_0)$ for each $w\in\mathcal{W}$.
We have
\begin{equation}\label{property1 for h}
\mathscr{D}_k h_{T,M}(\mu)\ll_k M^{-k},
\end{equation}
for any differential operator $\mathscr{D}_k$ of order $k$,
which we use frequently when we integrate by parts,
and sufficiently many differentiations can save arbitrarily many powers of $T$.
Moreover, by trivial estimate, we have
\begin{equation}\label{property2 for h}
  \int_{\Re \mu=0}h_{T,M}(\mu)\dd_{\spec}\mu\asymp T^3M^2.
\end{equation}

\subsection{The weight functions}

For the weight functions, we will need the following results
in Blomer--Buttcane \cite[Lemma 1, Lemma 8, and Lemma 9]{blomer2015subconvexity}.

\begin{lemma}\label{truncate 1st}
  For some large enough constant $B>0$, we have
  \[
    \Phi_{w_4}(y) \ll |y|^{1/10}T^B,\quad
    \Phi_{w_5}(y) \ll |y|^{1/10}T^B,\quad
    \Phi_{w_6}(y_1,y_2) \ll |y_1y_2|^{3/5}T^B.
  \]
\end{lemma}

\begin{lemma}\label{truncate 2nd}
  \begin{itemize}
    \item [(i)]    If $0<|y|\leq T^{3-\ve}$, then for any constant $B>0$, we have
                    \[
                      \Phi_{w_4}(y) \ll_{\ve,B} T^{-B}.
                    \]
    \item [(ii)]   If $T^{3-\varepsilon}<y$, then
                    \[
                      y^{k}\Phi_{w_4}^{(k)}(y)\ll_{k,\varepsilon} T^{3+2\varepsilon}M^2(T+|y|^{1/3})^{k}.
                    \]

  \end{itemize}
\end{lemma}

\begin{lemma}\label{truncate 3rd}
  Let $\Upsilon:=\min\{|y_1|^{1/3}|y_2|^{1/6},|y_1|^{1/6}|y_2|^{1/3}\}$.
  If $\Upsilon\ll T^{1-\ve}$, then we have
  \[ \Phi_{w_6}(y_1,y_2) \ll T^{-B}. \]
\end{lemma}

There is a slight difference that we use \eqref{property1 for h} and \eqref{property2 for h}
instead of \cite[(3.7) and (3.8)]{blomer2015subconvexity}.
This have no influence in the proof.
We define $\Phi_{2,w_4}(y)$, $\Phi_{2,w_5}(y)$, and $\Phi_{2,w_6}(y)$ as in \eqref{defPhi} by using
the test function $h_2(\mu)$.
In the proof of the above three Lemmas, the only two properties of $h_{T,M}(\mu)$ which are used is \eqref{property1 for h} and \eqref{property2 for h}.
Here, we remark that $h_2(\mu)$ also satisfies these two inequalities
by using properties of $W_{\mu,N}$ (see \S \ref{sec:2ndmoment}).
So, $\Phi_{2,w_4}(y)$, $\Phi_{2,w_5}(y)$, and $\Phi_{2,w_6}(y)$ also satisfies the corresponding bounds in the above three Lemmas.

\section{The mollified second moment}\label{sec:2ndmoment}
Let $G(s)=(\cos\f{\pi s}{A})^{-100A}$, where $A$ is a positive integer.
For $|L(\thf,\phi_j)|^2$, we will use the
following approximate functional equation.

\begin{lemma}\label{lemma:AFE2}
  Let $\phi_j$ be a $GL(3)$ Maass form with Langlands parameters $\mu_j \in \Lambda_{1/2}'$.
  We have
  \[
    |L(\thf,\phi_j)|^2 = 2 \sum_{m_1,m_2}
    \frac{\overline{A_j(m_1,m_2)}}{(m_1m_2)^{1/2}}W_{j}\left(m_1m_2\right),
  \]
  where
  \[
    W_{j}(y):= \frac{1}{2\pi i}\int_{(3)} \zeta(1+2s) (\pi^3y)^{-s} \prod_{k=1}^{3}\prod_\pm \frac{\Gamma\left(\frac{s+1/2\pm\mu_{j,k}}{2}\right)}{\Gamma\left(\frac{1/2\pm\mu_{j,k}}{2}\right)}
    G(s)\frac{\dd s}{s}.
  \]
  Moreover, we have
  \[
    y^{i}W_{j}^{(i)}(y) \ll_{i,B} \bigg(1+\frac{y}{\prod_{k=1}^{3}(1+|\mu_{j,k}|)}\bigg)^B,
  \]
  for any non-negative integer $i$, and any large positive integer $B$.
\end{lemma}

\begin{proof}
  See e.g. Iwaniec--Kowalski \cite[\S5.2]{iwaniec2004analytic}.
\end{proof}

Note that the sum in Proposition \ref{second moment} is essentially supported on the generic forms
which satisfy
\[
  |\mu_{j,k}|\asymp|\nu_{j,k}|\asymp T, \quad 1\leq k \leq 3.
\]
So we assume $\phi_j$ also satisfy the above relation.
By Stirling's formula, if $|s|\ll |z|^{1/2}$, then
$$
\f{\Gamma(z+s)}{\Gamma(z)}=z^s\left(1+\sum_{n=1}^N\f{P_n(s)}{z^{n}}+O\Big(\f{(1+|s|)^{2N+2}}{|z|^{N+1}}\Big)\right),
$$
for certain polynomials $P_n(s)$ of degree $2n$. Since $G(s)$ has exponential decay, we may truncate at
$|\Im s|\ll T^\varepsilon$ with only a small error in $W_{j}(y)$.
Based on the above arguments, together with $\mu_j \in \Lambda_{1/2}'$, we have
\begin{align*}
    & W_{j}(y) = \frac{1}{2\pi i} \int_{(3)} \zeta(1+2s) (\pi^3y)^{-s} \prod_{k=1}^{3} \prod_\pm \Big(\frac{1/2\pm\mu_{j,k}}{2}\Big)^{s/2} \\
    & \hskip 80pt \cdot \left(1+\sum_{n=1}^N\f{2^nP_n(s)}{(1/2\pm\mu_{j,k})^{n}}
    +O\Big(\f{(1+|s|)^{2N+2}}{|\mu_k|^{N+1}}\Big)\right)
    G(s)\frac{\dd s}{s}.
\end{align*}
On the other hand, by Hecke multiplicativity relations we have
\[
  \begin{split}
    |M_j|^2 & = \sum_{\ell_1}\sum_{\ell_2} \frac{1}{\ell_1^{1/2}\ell_2^{1/2}}
      x_{\ell_1}\overline{x_{\ell_2}} A_j(\ell_1,1)\overline{A_j(\ell_2,1)} \\
    & =\sum_{d}\sum_{\ell_1}\sum_{\ell_2} \frac{1}{d\ell_1^{1/2}\ell_2^{1/2}}
        x_{d\ell_1}\overline{x_{d\ell_2}} A_j(\ell_1,\ell_2),
  \end{split}
\]
By Lemma \ref{lemma:AFE2} and Hecke relations again, we have
\[
  \begin{split}
    & \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} |L(\thf,\phi_j)M_j|^2 \\
    & \hskip 30pt = 2 \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} \sum_{d}\sum_{\ell_1}\sum_{\ell_2} \frac{x_{d\ell_1}\overline{x_{d\ell_2}}}{d\ell_1^{1/2}\ell_2^{1/2}}
         A_j(\ell_1,\ell_2)
        \sum_{m_1,m_2}\frac{1}{m_1^{1/2}m_2^{1/2}}
        W_{j}\left(m_1m_2\right)\overline{A_j(m_1,m_2)} \\
    & \hskip 30pt = 2 \sum_{d}\frac{1}{d}\sum_{\ell_1,\ell_2}
        \frac{x_{d\ell_1}\overline{x_{d\ell_2}}}{\ell_1^{1/2}\ell_2^{1/2}}
        \sum_{m_1,m_2}\frac{1}{m_1^{1/2}m_2^{1/2}}
        \sum_{j}\frac{h_{T,M}(\mu_j) W_{j}\left(m_1m_2\right)}{\cN_j}
        A_j(\ell_1,\ell_2) \overline{A_j(m_1,m_2)}.
  \end{split}
\]
Let
$$
  W_{j,N}(y) = \frac{1}{2\pi i} \int_{(3)} \zeta(1+2s) (\pi^3y)^{-s}\prod_{k=1}^{3}\prod_\pm \left(\frac{1/2\pm\mu_{j,k}}{2}\right)^{s/2}
  \left(1+\sum_{n=1}^N\f{2^nP_n(s)}{(1/2\pm\mu_{j,k})^{n}}\right)
    G(s)\frac{\dd s}{s}.
$$
Then we can replace $W_{j}$ by $W_{j,N}$ with a negligible error if we choose $N$ to be large enough.
Now we use the Kuznetsov trace formula (see Lemma \ref{lemma: KTF}) with the test function
\begin{align}\label{test h2}
h_2(\mu)=h_{T,M}(\mu)W_{\mu,N}(m_1m_2),
\end{align}
where $W_{\mu,N}$ is defined the same as $W_{j,N}$ with $\mu$ replacing $\mu_j$.
It turns out we are led to estimate
\begin{equation}\label{use Ku 2}
  \begin{split}
      2 \sum_{d}\frac{1}{d}\sum_{\ell_1,\ell_2}
        \frac{x_{d\ell_1}\overline{x_{d\ell_2}}}{\ell_1^{1/2}\ell_2^{1/2}}
        \sum_{m_1,m_2}\frac{1}{m_1^{1/2}m_2^{1/2}}
        \Big(\Delta^{(2)}+\Sigma_4^{(2)}+\Sigma_5^{(2)}+\Sigma_6^{(2)}
        -\mathcal{E}_{\mathrm{max}}^{(2)}-\mathcal{E}_{\mathrm{min}}^{(2)}\Big),
  \end{split}
\end{equation}
where
\begin{displaymath}
    \begin{split}
      \Delta^{(2)} & := \delta_{\ell_1, m_1} \delta_{\ell_2, m_2}
            \frac{1}{192\pi^5} \int_{\Re \mu = 0} h_2(\mu) \dd_{\spec}\mu,\\
      \Sigma_{4}^{(2)}& := \sum_{\epsilon  = \pm 1} \sum_{\substack{D_2 \mid D_1\\  m_2 D_1= \ell_1 D_2^2}}
            \frac{\tilde{S}(-\epsilon \ell_2, m_2, m_1; D_2, D_1)}{D_1D_2}
            \Phi_{2,w_4}\!\left(  \frac{\epsilon m_1m_2\ell_2}{D_1 D_2} \right),  \\
      \Sigma_{5}^{(2)} & := \sum_{\epsilon  = \pm 1} \sum_{\substack{D_1 \mid D_2\\ m_1 D_2 = \ell_2 D_1^2}}
            \frac{\tilde{S}(\epsilon \ell_1, m_1, m_2; D_1, D_2) }{D_1D_2}
            \Phi_{2,w_5}\!\left( \frac{\epsilon \ell_1m_1m_2}{D_1 D_2}\right),\\
      \Sigma_6^{(2)} & := \sum_{\epsilon_1, \epsilon_2 = \pm 1} \sum_{D_1,  D_2  }
            \frac{S(\epsilon_2 \ell_2, \epsilon_1 \ell_1, m_1, m_2; D_1, D_2)}{D_1D_2}
            \Phi_{2,w_6}\!  \left( - \frac{\epsilon_2 m_1\ell_2D_2}{D_1^2},
            - \frac{\epsilon_1 m_2\ell_1D_1}{ D_2^2}\right),
    \end{split}
\end{displaymath}
with $\Phi_{2,w_4}(y)$, $\Phi_{2,w_5}(y)$, and $\Phi_{2,w_6}(y)$ defined as in \eqref{defPhi} by using
the new test function $h_2(\mu)$ given by \eqref{test h2}, respectively;
and
\begin{equation*}
    \begin{split}
       \cE_{\max}^{(2)} & :=  \sum_{g} \frac{1}{2\pi i} \int\limits_{\Re(\mu)=0}
            \frac{h_2(\mu+\mu_g,\mu-\mu_g,-2\mu)}{\cN_{\mu,g}}
             \overline{B_{\mu,g}(m_1, m_2)} B_{\mu,g}(\ell_1, \ell_2) \dd\mu, \\
       \cE_{\min}^{(2)} & := \frac{1}{24(2\pi i)^2} \iint\limits_{\Re(\mu)=0} \frac{h_2(\mu)}{\cN_{\mu}}
                        \overline{A_{\mu}(m_1, m_2)} A_{\mu}(\ell_1, \ell_2) \dd\mu_1 \dd\mu_2.
    \end{split}
\end{equation*}

\subsection{The diagonal term}

Note that we have $\ell_1=m_1$, and $\ell_2=m_2$.
Thus we infer that the diagonal term in \eqref{use Ku 2} is
\begin{align}\label{diagonal}
  \f{1}{96\pi^5}\sum_{d}\sum_{\ell_1}\sum_{\ell_2}\f{x_{d\ell_1}x_{d\ell_2}}{\ell_1\ell_2}
  \int_{\Re \mu = 0} h_{T,M}(\mu) W_{\mu,N}\left(\ell_1\ell_2\right) \dd_{\spec}\mu.
\end{align}
By the definition of $W_{\mu,N}$, we only need to deal with the leading term,
which contributes
\begin{align*}
  \f{1}{96\pi^5}\int_{\Re \mu = 0} h_{T,M}(\mu) S(\mu) \dd_{\spec}\mu,
\end{align*}
where
\begin{align*}
  S(\mu) :=\sum_{d} \f{\mu^2(d)}{d} \sum_{(\ell_1,d)=1}\sum_{(\ell_2,d)=1} \f{\mu(\ell_1)a_{d\ell_1}\mu(\ell_2)a_{d\ell_2}}{\ell_1\ell_2}
  \mathcal{W}_{\mu}(\ell_1\ell_2),
\end{align*}
with
\begin{align*}
  a_{\ell} & := \f{1}{2\pi i}\int_{(3)}\f{(L/\ell)^s}{s^2}\frac{\dd s}{\log L}
  =  \begin{cases}
        \f{\log(L/\ell)}{\log L},\ \ &\ell\leq L,
        \\\noalign{\vskip 1mm}
        0,\ \ \ &\ell>L,
     \end{cases} \\
  \mathcal{W}_{\mu}(y) & := \f{1}{2\pi i} \int_{(3)} \zeta(1+2s)\left(\f{\pi^3 y}{|\mu_1\mu_2\mu_3|}\right)^{-s}G(s)\f{\dd s}{s}.
\end{align*}
By moving the line of integration in $\mathcal{W}_{\mu}$ to left,
we pass the double pole at $s=0$. Hence, by the residue theorem, we infer that
\begin{align*}
  \mathcal{W}_{\mu}(\ell_1\ell_2)=c_1\log|\mu_1\mu_2\mu_3|+c_2\log \ell_1\ell_2+c_3+O(T^{-B}),
\end{align*}
where $c_1$, $c_2$, and $c_3$ are constants. The main contribution from $\mathcal{W}_{\mu}(\ell_1\ell_2)$
comes from the first two terms, which can be treated similarly. For convenience, we only give the details
of the first term. Note that the goal is to show
\begin{align}\label{goal}
  \sum_{d}\f{\mu^2(d)}{d} \sum_{(\ell_1,d)=1} \sum_{(\ell_2,d)=1} \f{\mu(\ell_1)a_{d\ell_1}\mu(\ell_2)a_{d\ell_2}}{\ell_1\ell_2}
  \int_{\Re \mu=0}h_{T,M}(\mu) \log|\mu_1\mu_2\mu_3|\dd_{\spec}\mu
  \ll T^{3}M^2.
\end{align}
For the $\ell_1$-sum, we have
\begin{align}\label{l1-sum}
  \sum_{(\ell_1,d)=1}\f{\mu(\ell_1)a_{d\ell_1}}{\ell_1}
  =\frac{1}{2\pi i} \int_{(2)}\prod_{p|d}\left(1-\f{1}{p^{1+s}}\right)^{-1}
  \f{(L/(d\ell_1))^s}{s^2\zeta(1+s)\log L}\dd s.
\end{align}
Recall that in the region $\Re s\geq1-c/\log(|\Im s|+3)$ (here $c$ is some positive constant),
$\zeta(s)$ is analytic except for a single pole at $s=1$, and has no zeros and satisfies
$\zeta(s)^{-1}\ll \log(|\Im s|+3)$, $\zeta'(s)/\zeta(s)\ll \log(|\Im s|+3)$
(see e.g. \cite[(3.11.7) and (3.11.8)]{titchmarsh1986theory}).
We move the line of integration in \eqref{l1-sum} to
\begin{equation}\label{defGamma}
  \mathcal{C}_{\varepsilon}
  := \big\{{\rm i}x : |x|\geq\varepsilon\big\}\cup
  \big\{\varepsilon \text{e}^{{\rm i}\vartheta} : \tfrac{\pi}{2}\leq\vartheta\leq\tfrac{3\pi}{2}\big\},
\end{equation}
and $\varepsilon$ is sufficiently small.
It follows that
\begin{align}\label{l1-moving}
  \eqref{l1-sum} = \f{1}{\log L} \prod_{p\mid d} \frac{p}{p-1}
  +\int_{\mathcal{C}_{\varepsilon}}\prod_{p|d}\left(1-\f{1}{p^{1+s}}\right)^{-1}
  \f{(L/(d\ell_1))^s}{s^2\zeta(1+s)\log L}\dd s.
\end{align}
We have the similar expression for the $\ell_2$-sum.
Inserting these into \eqref{goal}, we consider the resulting $d$-sum.
A typical term is
\begin{align}\label{typical term}
\sum_{d\leq L}\f{\mu^2(d)}{d}\prod_{p\mid d} \left(\frac{p}{p-1}\right)^2\ll \log T
\end{align}
which implies \eqref{goal} by trivial computation.
For the term involving $\log(\ell_1\ell_2)=\log \ell_1+\log \ell_2$,
we have
\begin{align*}
  \sum_{\substack{\ell_1\geq1\\(\ell_1,d)=1}}\f{\mu(\ell_1)\log \ell_1}{\ell_1^{1+s}}
  = - \bigg\{\f{1}{\zeta(1+s)}\prod_{p\mid d} \bigg(1-\frac{1}{p^{s+1}}\bigg)^{-1}\bigg\}'.
\end{align*}
A similar argument shows that its contribution to \eqref{diagonal} is $\ll T^3M^2$.

\subsection{The $w_4$ and $w_5$ terms}

We only deal with the $w_4$-term, since the $w_5$-term is very similar.
Inserting a smooth unity to $m_1$, $m_2$ sums, we are led to estimate
\begin{align}\label{w4 contribution}
\sum_{\epsilon=\pm1}\sum_{d}\frac{1}{d}\sum_{\ell_1,\ell_2}
        \frac{1}{\sqrt{\ell_1\ell_2}}
        \Sigma_4^{(2)}(\epsilon,d,\ell_1,\ell_2,M_1,M_2),
\end{align}
where $M_1,M_2\gg1$, $M_1M_2\ll T^{3+\varepsilon}$,
\begin{equation*}
  \begin{split}
    & \Sigma_4^{(2)}(\epsilon,d,\ell_1,\ell_2,M_1,M_2) \\
    & \hskip 30pt :=\sum_{m_1=1}^\infty\sum_{m_2=1}^\infty \f{W_1(\f{m_1}{M_1})W_2(\f{m_2}{M_2})}{\sqrt{m_1m_2}}
      \sum_{\substack{D,\delta\\ m_2\d=\ell_1D}}\frac{\tilde{S}(-\epsilon \ell_2, m_2, m_1; D, D\d)}{D^2\d}
            \Phi_{2,w_4}\!\left(  \frac{\epsilon m_1m_2\ell_2}{D^2\d} \right),
  \end{split}
\end{equation*}
and $W_i(x)$ ($i=1,2$) are compactly supported in $[1,2]$ and satisfy $x^jW_i^{(j)}(x)\ll 1$.
By Lemma \ref{truncate 1st}, we can truncate the $D$, $\d$ sums
at some $T^B$ for some larege $B$ at the cost of a negligible error.
Then by Lemma \ref{truncate 2nd}, we can truncate the sums, again with a negligible error, at
\begin{align}\label{eqn:condition1}
  \f{m_1m_2\ell_2}{D^2\d}\geq T^{3-\varepsilon},
\end{align}
or in other words,
$$
D^2\d\leq \f{M_1M_2L}{T^{3-\varepsilon}}\ll LT^\varepsilon.
$$
Note that we have $m_2\d=\ell_1D$ now, which implies $M_2$ is small. That is,
$$
  M_2 \leq LD/\d \leq (M_1M_2/T^{3-\varepsilon})^{1/2}L^{3/2}/\d^{3/2}
  \ll L^{3/2}T^\varepsilon/\d^{3/2}.
$$
And by \eqref{eqn:condition1} we have
$$
  M_1\gg T^{3-\varepsilon}D^2\delta/(m_2\ell_2)
  = T^{3-\varepsilon}D\delta^2/(\ell_1\ell_2)\gg T^{3-\varepsilon}/L^2.
$$
We apply the Poisson summation formula in the $m_1$-variable, getting
\begin{align*}
  & \sum_{m_1=1}^\infty \f{W_1(\f{m_1}{M_1})}{\sqrt{m_1}}\tilde{S}(-\epsilon \ell_2, m_2, m_1; D, D\d)
  \Phi_{2,w_4}\!\left(  \frac{\epsilon m_1m_2\ell_2}{D^2\d} \right)
  \\
  & \hskip 50pt =\sum_{a(\mod \d)}\tilde{S}(-\epsilon \ell_2, m_2, a; D, D\d)\sum_{m_1\equiv a(\mod \d)}
  \f{W_1(\f{m_1}{M_1})}{\sqrt{m_1}}\Phi_{2,w_4}\!\left(  \frac{\epsilon m_1m_2\ell_2}{D^2\d} \right)
  \\
  & \hskip 50pt =\f{M_1^{1/2}}{\d}\sum_{a(\mod \d)}\tilde{S}(-\epsilon \ell_2, m_2, a; D, D\d)
  \sum_{m\in\mathbb{Z}}e\left(\f{am}{\d}\right)\\
  & \hskip 150pt  \cdot
  \int_{-\infty}^\infty \f{W_1(x)}{x^{1/2}}
  \Phi_{2,w_4}\!\left(\frac{\epsilon m_2\ell_2M_1x}{D^2\d} \right)e\left(-\f{mM_1x}{\d}\right)\dd x.
\end{align*}
Integration by parts in connection with Lemma \ref{truncate 2nd} ${\rm (ii)}$ and the above bounds for $M_1$ and $M_2$
shows that the integral is negligible unless $m=0$.
When $m=0$, by opening $\tilde{S}(-\epsilon \ell_2, m_2, a; D, D\d)$ and compute the $a$ sum, we obtain
\begin{align*}
\sum_{a(\mod \d)}\tilde{S}(-\epsilon \ell_2, m_2, a; D, D\d)=\sum_{a(\mod \d)}\sum_{\substack {C_1(\mod D), C_2(\mod D\d)\\(C_1,D)=(C_2,\d)=1}}e\left(m_2\f{\bar{C_1}C_2}{D}+a\f{\bar{C_2}}{\d}-\epsilon \ell_2\f{C_1}{D}\right)
=0,
\end{align*}
unless $\d=1$.
With the help of this and $m_2=\ell_1D$, we see that the contribution from $m=0$ to
$\Sigma_4^{(2)}(\epsilon,d,\ell_1,\ell_2,M_1,M_2)$ is
\begin{align*}
M_1^{1/2}\ell_1^{-1/2}\sum_{D=1}^\infty W_2\left(\f{\ell_1D}{M_2}\right)\f{1}{D^{3/2}}\sum_{\substack{C_1(\mod D)\\(C_1,D)=1}}
e\left(-\f{\epsilon \ell_2C_1}{D}\right)
\int_{-\infty}^\infty \f{W_1(x)}{x^{1/2}}
\Phi_{2,w_4}\!\left(\frac{\epsilon \ell_1\ell_2M_1x}{D} \right)\dd x.
\end{align*}
By inserting the definitions of $\Phi_{2,w_4}(x)$ and $K_{w_4}(y;\mu)$, the above $x$-integral becomes
\begin{align}\label{inserting Phi}
  &\int_{-\infty}^\infty \f{W_1(x)}{x^{1/2}}\int_{\Re \mu=0}h_2(\mu)K_{w_4}\left(\frac{\epsilon \ell_1\ell_2M_1x}{D};\mu\right)\dd_{\spec}\mu \dd x\nonumber
\\
  &\hskip 30pt =\int_{\Re \mu=0}h_2(\mu)\int_{-i\infty}^{i\infty}\left(\int_{-\infty}^\infty x^{-s-1/2}W_1(x)\dd x\right)
\left|\frac{\ell_1\ell_2M_1}{D}\right|^{-s}\tilde{G}^\epsilon(s,\mu)\f{ds}{2\pi i}\dd_{\spec}\mu\nonumber
\\
  &\hskip 30pt =\int_{\Re \mu=0}h_2(\mu)\int_{-i\infty}^{i\infty}\left|\frac{\ell_1\ell_2M_1}{D}\right|^{-s}
\hat{W}_1\left(\f{1}{2}-s\right)\tilde{G}^\epsilon(s,\mu)\f{\dd s}{2\pi i}\dd_{\spec}\mu,
\end{align}
where the Mellin transform $\hat{W}_1$ is entire and rapidly decaying, which means that
we can restrict the $s$-integral to $|\Im s|\leq T^\varepsilon$.
Inserting the definition of $\tilde{G}^\epsilon(s,\mu)$ into \eqref{inserting Phi}, we obtain
the corresponding $\mu$-integral:
\begin{align*}
\int_{\Re \mu=0}h(\mu)\Biggl(\prod_{j=1}^3\frac{\Gamma(\frac{1}{2}(s-\mu_j))}{\Gamma(\frac{1}{2}(1-s+\mu_j))} \pm i   \prod_{j=1}^3\frac{\Gamma(\frac{1}{2}(1+s-\mu_j))}{\Gamma(\frac{1}{2}(2-s+\mu_j))} \Biggr)\dd_{\spec}\mu.
\end{align*}
We only need to consider the first part, since the second part is very similar.
For fixed $\sigma\in\mathbb{R}$, $t\in\mathbb{R}$ and $|t|\geq 1$, we have the Stirling formula
\begin{align*}
  \Gamma(\sigma+it)=\sqrt{2\pi}e^{\pi|t|/2}|t|^{\sigma-1/2}e^{it(\log|t|-1)} e^{i(\sigma-1/2)\lambda\pi/2} \left(1+O(|t|^{-1})\right),
\end{align*}
where $\lambda=1$, if $t\geq1$, and $\lambda=-1$, if $t\leq -1$.
For convenience, we denote $\mu_1=it_1$, $\mu_2=it_2$, $\mu_3=-i(t_1+t_2)=it_3$, and $\Im s=t$.
Here, in the essential integrated range, we have $|t_k|\asymp|t_k-t_k'|\asymp T$ for $1\leq k, k'\leq 3$ and $k\neq k'$,
since we are considering the generic case.
Without loss of generality, we assume $t_1>0$, $t_2>0$, and $t_3<0$, and get
\begin{equation}\label{use stirling}
  \begin{split}
     & \quad\ \prod_{j=1}^3\frac{\Gamma(\frac{1}{2}(s-\mu_j))}{\Gamma(\frac{1}{2}(1-s+\mu_j))} \\
     & = 2^{-3/2}e^{i(-3t(\log 2+1)+\pi/2)}
       \f{e^{i(t-t_1)\log(t_1-t)+i(t-t_2)\log(t_2-t)+i(t_1+t_2+t)\log(t_1+t_2+t)}}
       {(t_1-t)^{1/2}(t_2-t)^{1/2}(t_1+t_2+t)^{1/2}}
       + O(T^{-5/2}).
  \end{split}
\end{equation}
By trivial estimate, the error term from the above to \eqref{w4 contribution}
is $O(T^{2+\varepsilon}M^2L)$.

For the main term in \eqref{use stirling}, we use partial integration to prove
its contribution is small. Actually, for the $\mu_1$-integral, we denote
the phase function
$$
\phi(t_1):=(t-t_1)\log(t_1-t)+(t_1+t_2+t)\log (t_1+t_2+t).
$$
We see that
$$
\phi'(t_1)=\log\f{t_1+t_2+t}{t_1-t}\gg 1.
$$
And hence, by partial integration many times, the integral is negligible.
We finally prove that the contribution from the $w_4$-term is $O(T^{2+\varepsilon}M^2L)$.

\subsection{The $w_6$ term}

By Lemma \ref{truncate 1st}, we can truncate the $D_1$, $D_2$ sums at some $T^B$ for some large $B$ at the cost of a negligible error.
Then by Lemma \ref{truncate 3rd}, we can truncate the sum further at
\[
  \f{(m_1\ell_2)^{1/3}(m_2\ell_1)^{1/6}}{D_1^{1/2}}\geq T^{1-\varepsilon},
  \quad \textrm{and} \quad
  \f{(m_1\ell_2)^{1/6}(m_2\ell_1)^{1/3}}{D_2^{1/2}}\geq T^{1-\varepsilon},
\]
which means
$$
\f{(m_1m_2\ell_1\ell_2)^{1/2}}{(D_1D_2)^{1/2}}\geq T^{2-\varepsilon}.
$$
This gives us
$$
D_1D_2\ll \f{m_1m_2\ell_1\ell_2}{T^{4-\varepsilon}},
$$
which is impossible provided that $L\leq T^{1/2-\varepsilon}$, since the essential sums of $m_1$ and $m_2$
are truncated at $m_1m_2\leq T^{3+\varepsilon}$.
Thus, the contribution of the $w_6$ term is $O(T^{-B})$.

\subsection{The contribution from the Eisentein series}

We only treat the contribution of the maximal Eisenstein series,
since the minimal Eisenstein series can
be handled similarly and the contribution will be smaller.
We have the Weyl law on $GL_2$ (see \cite{hejhal}),
$$
\sharp\left\{g: t_g\leq T\right\}=\frac{T^2}{12}-\frac{T\log T}{2\pi}+C_0T+O\left(\frac{T}{\log T}\right),
$$
where $C_0$ is a constant.
Combining this together with definitions of $h(\mu)$ and $h_2(\mu)$, and note that
we are considering the generic case, we see that the essential region of the integration
on $\mu$ is of length $\asymp M$, and the essential number of the sum of $\mu_g$ is of
size $\asymp TM$. Hence, by the bound \eqref{eqn:Amax<<}, \eqref{N max} and the above argument,
we have
\[
  \mathcal{E}_{\mathrm{max}}^{(2)} \ll_{\varepsilon}
  TM^2  (T^3L^2)^{7/64} T^{\varepsilon}.
\]
Hence the contribution of the right hand side of \eqref{use Ku 2}
of the maximal Eisenstein series is
\[
  \ll_{\varepsilon} L T^{3/2} TM^2  (T^3L^2)^{7/64} T^{\varepsilon}
  \ll_{\varepsilon} T^{3}M^2 T^{-11/64+\varepsilon}L^{39/32}
  \ll_{\varepsilon} T^{3}M^2,
\]
provided $L\ll T^{11/78-\varepsilon}$.

\section{The mollified first moment}\label{sec:1stmoment}

Let $G(s)$ be defined as in \S\ref{sec:2ndmoment}, and $T_0=T^{1+\varepsilon}$.
Consider the integral
\begin{align*}
  \f{1}{2\pi i}\int_{(3)}G(s)L\left(s+\f{1}{2},\phi_j\right)\f{T_0^{3s}}{s}\dd s.
\end{align*}
By moving the line of integration to $\Re(s)=-3$, and using the functional equation \eqref{eqn:FE}, we get
\begin{align*}
  L\left(s+\f{1}{2},\phi_j\right)=\sum_{m\geq 1}\f{A_j(1,m)}{m^{1/2}}V(m,T_0)
  + \sum_{m\geq 1}\f{A_j(m,1)}{m^{1/2}}V_{j}(m,\mu_j),
\end{align*}
where
$$
  V(y,T_0):=\f{1}{2\pi i}\int_{(3)}G(s)\left(\f{y}{T_0^3}\right)^{-s}\f{\dd s}{s},
$$
and
$$
  V_{j}(y,\mu_j):=\int_{(3)}G(s)(yT_0^3)^{-s} \prod_{k=1}^3\f{\Gamma\left(\f{s+1/2+\mu_{j,k}}{2}\right)}{\Gamma\left(\f{-s+1/2-\mu_{j,k}}{2}\right)} \f{\dd s}{s}.
$$
To see the properties of $V$ and $V_{j}$, we use the strategy in \cite{iwaniec2004analytic}.
Obviously, we have
$V(y,T_0)=1+O\left((y/T_0^3)^B\right)$,
and
$V(y,T_0)\ll (y/T_0^3)^{-B}$,
by moving the integration line to $\Re s=-B$ and $\Re s=B$ respectively.
For $V_{j}$, we have $V_{j}(y,\mu_j)\ll (yT)^{-B}$ for any $y\geq 1$,
by using the Stirling formula and moving the integration line
to $\Re s=B/(3\varepsilon)$.
With the help of these, we infer that
\[\label{first moment}
  \begin{split}
    & \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} L(\thf,\phi_j) M_j \\
    & \hskip 50pt = \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} \sum_{\ell} \frac{x_{\ell}}{\ell^{1/2}}
         A_j(1,\ell)
        \sum_{m}\frac{1}{m^{1/2}}
        V\left(m,T_0\right)\overline{A_j(m,1)} + O(T^{-B})\\
    & \hskip 50pt =  \sum_{\ell}
        \frac{x_{\ell}}{\ell^{1/2}}
        \sum_{m}\frac{V(m,T_0)}{m^{1/2}}
        \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j}
        A_j(1,\ell) \overline{A_j(m,1)}+O(T^{-B}).
  \end{split}
\]
Then, by the Kuznetsov trace formula, we have
\[
  \begin{split}
   & \sum_{\ell}\frac{x_{\ell}}{\ell^{1/2}}
     \sum_{m}\frac{V(m,T_0)}{m^{1/2}}
     \sum_{j} \frac{h_{T,M}(\mu_j)}{\cN_j} A_j(1,\ell) \overline{A_j(m,1)}\\
   & \hskip 50pt = \sum_{\ell}\frac{x_{\ell}}{\ell^{1/2}}
     \sum_{m}\frac{V(m,T_0)}{m^{1/2}}
     \left(\Delta^{(1)}+\Sigma_4^{(1)}+\Sigma_5^{(1)}+\Sigma_6^{(1)}
     -\mathcal{E}_{\mathrm{max}}^{(1)}-\mathcal{E}_{\mathrm{min}}^{(1)}\right),
  \end{split}
\]
where
\begin{displaymath}
    \begin{split}
      \Delta^{(1)} & := \delta_{\ell, 1} \delta_{1, m}
            \frac{1}{192\pi^5} \int_{\Re \mu = 0} h_{T,M}(\mu) \dd_{\spec}\mu,\\
      \Sigma_{4}^{(1)}& := \sum_{\epsilon  = \pm 1} \sum_{\substack{D_2 \mid D_1\\  m D_1= n_1 D_2^2}}
            \frac{\tilde{S}(-\epsilon, m, 1; D_2, D_1)}{D_1D_2}
            \Phi_{w_4}\!\left(  \frac{\epsilon m}{D_1 D_2} \right),  \\
      \Sigma_{5}^{(1)} & := \sum_{\epsilon  = \pm 1} \sum_{\substack{D_1 \mid D_2\\ D_2 = D_1^2}}
            \frac{\tilde{S}(\epsilon \ell, 1, m; D_1, D_2) }{D_1D_2}
            \Phi_{w_5}\!\left( \frac{\epsilon \ell m}{D_1 D_2}\right),\\
      \Sigma_6^{(1)} & := \sum_{\epsilon_1, \epsilon_2 = \pm 1} \sum_{D_1,  D_2  }
            \frac{S(\epsilon_2, \epsilon_1 \ell, 1, m; D_1, D_2)}{D_1D_2}
            \Phi_{w_6}\!  \left( - \frac{\epsilon_2 D_2}{D_1^2},
            - \frac{\epsilon_1 m \ell D_1}{ D_2^2}\right),
    \end{split}
\end{displaymath}
and
\begin{equation*}
    \begin{split}
       \cE_{\max}^{(1)} & :=  \sum_{g} \frac{1}{2\pi i} \int\limits_{\Re(\mu)=0}
            \frac{h_{T,M}(\mu+\mu_g,\mu-\mu_g,-2\mu)}{\cN_{\mu,g}}
             \overline{B_{\mu,g}(1, m)} B_{\mu,g}(\ell,1) \dd\mu, \\
       \cE_{\min}^{(1)} & := \frac{1}{24(2\pi i)^2} \iint\limits_{\Re(\mu)=0} \frac{h_{T,M}(\mu)}{\cN_{\mu}}
                        \overline{A_{\mu}(1, m)} A_{\mu}(\ell,1) \dd\mu_1 \dd\mu_2.
    \end{split}
\end{equation*}

\subsection{The diagonal term}

By trivial estimation, the contribution of the diagonal term is
\begin{align*}\label{diagonal first moment}
  V(1,T_0)\f{1}{192\pi^5}\int_{\Re \mu = 0} h_{T,M}(\mu) \dd_{\spec}\mu
  =\f{1}{192\pi^5}\int_{\Re \mu = 0} h_{T,M}(\mu) \dd_{\spec}\mu+O(T^{-B})\asymp T^3M^2.
\end{align*}

\subsection{The other terms}

The treatments of the other terms are actually similar and in fact easier than those
in the mollified second moment. So we omit the arguments here.

redpajama_set_name: RedPajamaArXiv