\section{Introduction}
Recently, there has been a surge in the interest of acquiring a theoretical understanding over deep neural network's behavior. Breakthroughs have been made in characterizing the optimization process, showing that learning algorithms such as stochastic gradient descent (SGD) tend to end up in one of the many local minima which have close-to-zero training loss \citep{choromanska2015loss, dauphin2014identifying, kawaguchi2016deep, nguyen2018optimization, du2018gradient}. However, these numerically similar local minima typically exhibit very different behaviors in terms of generalizability. It is, therefore, natural to ask two closely related questions: (a) What kind of local minima can generalize better? (b) How to find those better local minima?

To our knowledge, existing work focused only on one of the two questions. For the ``what'' question, various definitions of ``flatness/sharpness'' have been introduced and analyzed \citep{keskar2016large, neyshabur2017pac, neyshabur2017exploring, wu2017towards, liang2017fisher}. However, they suffer from one or more of the problems: (1) being mostly theoretical with no or poor empirical evaluations on modern neural networks, (2) lack of theoretical analysis and understanding, (3) in practice not applicable to finding better local minima. Regarding the ``how'' question, existing approaches ~\citep{hochreiter1997flat, sokolic2017robust, chaudhari2016entropy, hoffer2017train, neyshabur2015path, izmailov2018averaging} share some of the common drawbacks: (1) derived only from intuitions but no specific metrics provided to characterize local minima, (2) no or weak analysis of such metrics, (3) not applicable or no consistent generalization improvement for modern DNNs.

In this paper, we tackle both the ``what'' and the ``how'' questions in a unified manner. Our answer provides both the theory and applications for the generalization problems across different local minima. Based on the determinant of Fisher information estimated from the training set, we propose a metric that \emph{solves all the aforementioned issues}. The metric can well capture properties that characterize local minima of different generalization ability. We provide its theoretical analysis, primarily a generalization bound based on PAC-Bayes \citep{mcallester1999some, mcallester1999pac}. For modern DNNs in practice, it is necessary to provide a tractable approximation of our metric. We propose an intuitive and efficient approximation to compare it across different local minima. Our empirical evaluations fully illustrate the effectiveness of the metric as a strong indicator of local minima's generalizability. Moreover, from the metric we further derive and  design a practical regularization technique that guides the optimization process in finding better generalizable local minima. The experiments on image classification datasets demonstrate that our approach gives consistent generalization boost for a range of DNN architectures. Codes are available at \url{https://github.com/SeanJia/InfoMCR}.

\section{Related Work} \label{related_work}

It has been empirically shown that larger batch sizes lead to worse generalization \citep{keskar2016large}. \citet{hoffer2017train} analyzed how the training dynamics is affected by different batch sizes and presented a perturbed batch normalization technique for better generalization. While it effectively improves generalization for large-batch training, a specific metric that indicates the generalizability is missing. Similarly, \citet{elsayed2018large} employed a structured margin loss to improve performance of DNNs w.r.t. noise and adversarial attack yet no metric was proposed. Furthermore, this approach essentially provided no generalization gain in the normal training setup.

The local entropy of the loss landscape was proposed to measure ``flatness'' in \citet{chaudhari2016entropy}, which also designed an entropy-guided SGD that achieves faster convergence in training DNNs. However, the method does not consistently improve generalization, e.g., a decrease of performance on CIFAR-10 \citep{krizhevsky2009learning}. Another method that focused on modifying the optimization process is the Path-SGD proposed by \citet{neyshabur2015path}. Specifically, the authors derived an approximate steepest descent algorithm that utilizes the path-wise norm regularization to achieve better generalization. The authors only evaluated it on a two-layer neural network, very likely since the path norm is computationally expensive to optimize during training.

A flat minimum search algorithm was proposed by \citet{hochreiter1997flat} based on the ``flatness'' of local minima defined as the volume of local boxes. Yet since the boxes have their axes aligned to the axes of the model parameters, their volumes could be significant underestimations of ``flatness'' for over-parametrized networks, due to the specific spectral density of Hessian of DNNs studied in \citet{pennington2018spectrum, sagun2017empirical}. The authors of \citet{wu2017towards} also characterized the ``flatness'' by volumes. They considered the inverse volume of the basin of attraction and proposed to use the Frobenius norm of Hessian at the local minimum as a metric. In our experiments, we show that their metric does not accurately capture the generalization ability of local minima under different scenarios. Moreover, they have not derived a  regularizer from their metric.

Based on a ``robustness'' metric, \citet{sokolic2017robust} derived a regularization technique that successfully improves generalization on multiple image classification datasets. Nevertheless, we show that their metric fails to capture the generalizability across different local minima.

By using the Bayes factor, \citet{mackay1992practical} studied the generalization ability of different local minima obtained by varying the coefficient of L2 regularization. It derived a formula involving the determinant of Hessian, similar to the one in ours. Whereas, this approach has restricted settings and, without proposing an efficient approximation, its metric is not applicable to modern DNNs, let alone serving as a regularizer. A generalization bound is missing in \citet{mackay1992practical} as well.

In a broader context of the ``what'' question, properties that capture the generalization of neural networks have been extensively studied. Various complexity measures for DNNs have been proposed based on norm, margin, Lipschitz constant, compression and robustness \citep{bartlett2002rademacher, neyshabur2015norm, sokolic2017robust, xu2012robustness, bartlett2017spectrally, zhou2018nonvacuous, dziugaite2017computing, arora2018stronger, jiang2018predicting}. While some of them aimed to provide tight generalization bounds and some of them to provide better empirical results, none of the above approaches explored the ``how'' question at the same time.

Very recently, \citet{karakida2019universal} and \citet{sun2019lightlike} studied the Fisher information of the neural network through the lens of its spectral density. In specific, \citet{karakida2019universal} applied mean-field theory to study the statistics of the spectrum and the appropriate size of the learning rate. Also, an information-theoretic approach, \citet{sun2019lightlike} derived a novel formulation of the minimum description length in the context of deep learning by utilizing tools from singular semi-Riemannian geometry.

\section{Outline and Notations} \label{pform}
In a typical $K$-way classification setting, each sample $x \in \mathcal{X}$ belongs to a single class denoted $c_x \in \{0, 1, ..., K\}$ according to the probability vector $y \in \mathcal{Y}$, where $\mathcal{Y}$ is the k-dimensional probability simplex so that $p(c_x=i) = y_i$ and $\sum_i y_i = 1$. Denote a feed-forward DNN parametrized by $w \in \mathbb{R}^W$ as $f_w: ~\mathcal{X} \rightarrow \mathcal{Y}$, which uses nonlinear activation functions and a softmax layer at the end. Denote the cross entropy loss as $\ell(f_w(x), y) = - \sum_i y_i \ln f_w(x)_i$. Denote the training set as $\mathcal{S}$, defined over $\mathcal{X} \times \mathcal{Y}$ with $|\mathcal{S}| = N$. The training objective is given as $\mathcal{L}(\mathcal{S}, w) = \frac{1}{N}\sum_{(x, y) \sim \mathcal{S}} \ell(f_w(x), y)$. Assume $\mathcal{S}$ is sampled from some true data distribution denoted $\mathcal{D}$, we can define expected loss $\mathcal{L}(\mathcal{D}, w) = \mathop{\mathbb{E}}_{(x, y) \sim \mathcal{D}}[\ell(f_w(x), y)]$. Throughout this paper, we refer a local minimum of $\mathcal{L}(\mathcal{S}, w)$ corresponding to a local minimizer $w_0$ as just the local minimum $w_0$. Our paper's outline and main achievements are:
\begin{itemize}
    \item In Sec. \ref{fisher} we relates Fisher information to neural network training as a prerequisite.
    \item In Sec. \ref{metric} we propose a metric $\gamma(w_0)$ that well captures local minima's generalizability.
    \item In Sec. \ref{gbound} we provide a generalization bound related to $\gamma(w_0)$.
    \item In Sec. \ref{appro} we propose an approximation $\widehat{\gamma}(w_0)$ for $\gamma(w_0)$, which is shown to be very effective in Sec. \ref{exp_char} via extensive empirical evaluations.
    \item In Sec. \ref{regularization} we devise a practical regularizer from $\widehat{\gamma}(w_0)$ that consistently improves generalizability across different DNNs, as evaluated in Sec. \ref{exp_reg}.
\end{itemize}

\subsection{Other Notations} \label{nota}
Denote $\nabla_w$ as gradient, $\mathbf{J}_w[\cdot]$ as Jacobian matrix, $\nabla_w^2$ as Hessian, $D_{\mathrm{KL}}(\cdot \Vert \cdot)$ as KL divergence, $\norm{\cdot}_2$ as spectrum or Euclidean norm, $\norm{\cdot}_F$ as Frobenius norm, $|\cdot|$ as determinant, $\textrm{tr}(\cdot)$ as trace norm, $\rho(\cdot)$ as spectral radius, $\ell\ell_{\mathcal{S}}(w)$ as log-likelihood on $\mathcal{S}$, and $[\cdot]_i$ for selecting the $i^{\rm{th}}$ entry.

We define $\boldsymbol{\ell}_x(w) \in \mathbb{R}^K$ whose $i^{\rm{th}}$ entry is $-\ln f_w(x)_i$ so that $\ell(f_w(x), y) = \boldsymbol{\ell}_x(w)^T y$. We define $\tilde{y} \in \mathbb{R}^K$ as the one-hot version of $y$, i.e., only keep the largest dimension as 1. Then we define $\tilde{\mathcal{L}}(\mathcal{S}, w) \in \mathbb{R}^N$ as the one-hot and vectorized version of $\mathcal{L}(\mathcal{S}, w)$, i.e., a vector whose entries are $\ell(f_w(x), {\tilde{y}})$ for $(x, y) \in \mathcal{S}$.
In other words, we approximate the cross entropy loss $\ell(f_w(x), y)$ by $\ell(f_w(x), \tilde{y})$.

\section{Local Minimum and Fisher Information} \label{fisher}
First of all, if $y$ is strictly one-hot and the training accuracy achieved at $w_0$ is 100\%, then $w_0$ cannot be a local minimizer, because the cross entropy loss remains positive even if arbitrarily close to zero. To admit local minima of full training accuracy, we assume the widely used label smoothing (LS) \citep{szegedy2016rethinking} is applied to train all models in our analysis. LS enables us to assume a local minimum $w_0$ of the training loss with $\sum_{(x, y) \in \mathcal{S}} D_{\mathrm{KL}}(f_{w_0}(x) \Vert y) = 0$. Although empirically we find that both our proposed metric and derived regularizer work similarly well without LS.

With LS in mind, each sample $(x, y) \in \mathcal{S}$ has its label $c_x$ sampled by $p(c_x=~i|x) = y_i$, denoted as $c_x  \sim y$. We denote the training data distribution as $(x, c_x) \sim \mathcal{S}$. The joint probability $p(x, c_x)$ modeled by the DNN is $p(x, c_x=i; w) = p(c_x=~i | x; w)\ p(x) = [f_w(x)]_i\ p(x)$ with $p(x) = \frac{1}{N}$. We can relate the training loss $\mathcal{L}(\mathcal{S}, w)$ to the negative log-likelihood $- \ell\ell_{\mathcal{S}}(w)$ by:
\begin{align*} \label{eq:1}
\mathcal{L}(\mathcal{S}, w) &= \frac{1}{N} \sum_{(x, y) \in \mathcal{S}} \boldsymbol{\ell}_x(w)^T y \\
&= - \frac{1}{N}\sum_{(x, y) \in \mathcal{S}} \mathop{\mathbb{E}}_{c_x \sim y} \ln p(c_x|x; w) \\
&= - \frac{1}{N} \ell\ell_{\mathcal{S}}(w) + \ln\frac{1}{N}
\end{align*}
\vspace{-20px}
\[\textrm{where}\ - \ell\ell_{\mathcal{S}}(w) = - \sum_{(x, y) \in \mathcal{S}} \mathbb{E}_{c_x \sim y} \ln p(x, c_x; w)\]
Also, $w_0$ corresponds to a local maximum of the likelihood function. The observed Fisher information \citep{efron1978assessing} evaluated at $w_0$ is defined using the Hessian of the negative log-likelihood, i.e.,
\begin{align} \label{eq:b}
    \mathcal{I}_{\mathcal{S}}(w_0) &= - \frac{1}{N} \nabla_w^2 \ell\ell_{\mathcal{S}}(w_0) =  \nabla_w^2 \mathcal{L}(\mathcal{S}, w_0) \nonumber \\
    &= \mathlarger{\mathop{\mathbb{E}}}_{(x, c_x) \sim \mathcal{S}} [\nabla_w \ln p_{w_0}(c_x) \nabla_w \ln p_{w_0}(c_x)^T]
\end{align}
where $p_{w_0}(c_x)$ denotes $p(c_x|x;w_0)$. The first equality is straightforward; the second has its proof in Appendix A. Since $p(c_x=i|x) = y_i$ and $\ln p(c_x =~ i|x; w_0) = [\boldsymbol{\ell}_x(w_0)]_i$, we can further simplify the Equation \ref{eq:b} to:
\begin{equation} \label{eq:1}
\mathcal{I}_{\mathcal{S}}(w_0) = \frac{1}{N}\sum_{(x, y) \in \mathcal{S}} \sum_{i=1}^K \nabla_{w}[\boldsymbol{\ell}_x(w_0)]_i \nabla_{w}[\boldsymbol{\ell}_x(w_0)]_i^T
\end{equation}
\textbf{Remark}: A global minimum $w_0$, if exists, is equivalent to a local minimum with 100\% training accuracy. At such $w_0$, we have $\nabla_w \ell(f_{w_0}(x), y) = \boldsymbol{0}$ as $D_{\mathrm{KL}}(f_{w_0}(x) \Vert y) = 0$; however, we also have $\mathcal{I}_{\mathcal{S}}(w_0) \in \mathbb{R}^{W\times W} \neq \boldsymbol{0}$.

\section{Local Minima Characterization} \label{charactering}

In this section, we derive and propose our metric, provide a PAC-Bayes generalization bound, and lastly, propose and give intuitions of an effective approximation of our metric for modern DNNs.

\subsection{Fisher Determinant as Generalization Metric} \label{metric}
We would like a metric to compare different local minima.
Under the Assumption \ref{assume_unique}, we can partition the parameter space of the neural network $f_w$ into disjoint regions, each is a small neighborhood of a local minimum taken into account. Formally, for a local minimum $w_0$ and a sufficiently small $V > 0$, we define the model class $\mathcal{M}(w_0)$ as the largest connected subset of $\{w \in ~\mathbb{R}^W: \mathcal{L}(\mathcal{S}, w) \leq h \}$ that contains $w_0$, where the height $h$ is defined as a real number such that the volume (namely the Lebesgue measure) of $\mathcal{M}(w_0)$ is $V$. By the Intermediate Value Theorem, for any sufficiently small $V$ there exists a corresponding height $h$. In essence, a local minimum $w_0$ of the entire parameter space becomes the global minimum of the model class $\mathcal{M}(w_0)$.

Formulated as a model class selection problem, we can compare different local minima by comparing their associated model classes.
We propose our metric $\gamma(\cdot)$, where lower $\gamma(w_0)$ indicates a better generalizable local minimum $w_0$:
\begin{equation}
\gamma(w_0) = \ln |\mathcal{I}_{\mathcal{S}}(w_0)|
\end{equation}
As a metric, $\gamma(w_0)$ requires $|\mathcal{I}_{\mathcal{S}}(w_0)|\ \neq 0$. Therefore, we state the following Assumption \ref{assume_unique}.

\begin{assm} \label{assume_unique}
The local minima $w_0$ we care about in the comparison are well isolated and unique in their corresponding neighborhood $\mathcal{M}(w_0)$.
\end{assm}

The Assumption \ref{assume_unique} is quite reasonable. For state-of-the-art network architectures used in practice, this is often the fact. To be precise, the Assumption \ref{assume_unique} is violated when the Hessian matrix at a local minimum is singular. Specifically, \citet{orhan2018skip} summarizes three sources of the singularity: (i) due to a dead neuron, (ii) due to identical neurons, and (iii) linear dependence of the neurons. As well demonstrated in \citet{orhan2018skip}, network with skip connection, e.g. ResNet \citep{he2016deep}, WRN \citep{wideresnet}, and DenseNet \citep{huang2017densely} used in our experiments, can effectively eliminate all the aforementioned singularity.

In \citet{dinh2017sharp}, the authors pointed out another source of the singularity specifically for networks with scale-invariant activation functions, e.g. ReLU. Namely, one can rescale the model parameters layer-wise so that the underlying function represented by the network remains unchanged in the region. In practice, this issue is not critical. Firstly, most modern deep ReLU networks, e.g. ResNet, WRN, and DenseNet, have normalization layers, e.g. BatchNorm \citep{ioffe2015batch}, applied before the activations. BatchNorm shifts all the inputs to the ReLU function, equivalently shifting the ReLU horizontally which makes it no longer scale-invariant. Secondly, due to the ubiquitous use of Gaussian weights initialization scheme and weight decay, most local minima obtained by gradient learning have weights of a relatively small norm. Consequently, in practice, we will not compare two local minima essentially the same but have one as the rescaled version of the other with a much larger norm of the weights.

Note that normally we have a limited size of the dataset, and so an approximation of $\gamma(w_0)$ is a must. We present our approximation scheme and its intuition in Sec. \ref{appro}.

\subsubsection{Connection to Fisher Information Approximation (FIA) Criterion}

Our metric $\gamma(w_0)$ is closely related to the FIA criterion. Based on the MDL principle \citep{rissanen1978modeling}, \citet{rissanen1996fisher} derived the FIA criterion to compare statistical models. Tailored to our setting, each model class $\mathcal{M}(w_0)$ has its FIA criterion as (lower FIA is better):
\begin{align*}
    \textrm{FIA} = &-\sum_{(x, y) \in \mathcal{S}} \mathlarger{\mathop{\mathbb{E}}}_{c_x \sim y} \ln p(x, c_x; w_0) \\
     &+ \frac{W}{2}\ln\frac{N}{2\pi} + \ln \int_{\mathcal{M}(w_0)} \sqrt{|\mathcal{J}(w)|}\ dw
\end{align*}
Where $\mathcal{J}(w)$ is the expected Fisher information evaluated at $w$. Notice that all regularity conditions of the FIA criterion are satisfied for the local minimum $w_0$ (also the global optimum of the model class), provided 100\% training accuracy and the Assumption \ref{assume_unique}. Ignoring the constant terms and assuming the training loss is locally quadratic in $\mathcal{M}(w_0)$ (later formalized and validated as Assumption \ref{assume}), the RHS becomes $\ln V + \frac{1}{2}\ln |\mathcal{J}(w_0)|$. Remind that $V$ is defined as the volume of $\mathcal{M}(w_0)$, also a constant.

Essentially in our metric we use the observed Fisher information in place of the expected one, making our metric tractable and applicable to modern DNNs.

\subsubsection{ Connection to Existing Flatness/Sharpness Metrics}

As mentioned in Sec. \ref{related_work}, the ``flatness'' of a local minimum was firstly related to the generalization ability of the neural network in \citet{hochreiter1997flat}, where the concept and the method are both preliminary. The idea is recently popularized in the context of deep learning by a series of paper such as \citet{keskar2016large, chaudhari2016entropy, wu2017towards}. Our approach roughly shares the same intuition with these existing works, namely, a ``flat'' local minimum admits less complexity and so generalizes better than a ``sharp'' one. To our best knowledge, our paper is the first among these work that provides both the theoretical analysis including a generalization bound and the empirical verification of both an efficient metric and a practical regularizer for modern network architectures.

\subsection{Generalization Bound} \label{gbound}

\begin{assm} \label{assume}
Given the training loss $\mathcal{L}(\mathcal{S}, w)$, its local minimum $w_0$ satisfying Assumption \ref{assume_unique} and the associated neighborhood $\mathcal{M}(w_0)$ whose volume $V$ is sufficiently small, as described in Sec. \ref{pform}, \ref{fisher} and \ref{metric}, respectively, when confined to $\mathcal{M}(w_0)$, we assume that $\mathcal{L}(\mathcal{S}, w)$ is quadratic.
\end{assm}

The Assumption \ref{assume} is quite reasonable as well. \citet{grunwald2007minimum} suggests that, a log-likelihood function, under regularity conditions (1) existence of its $1^{\rm{st}}$, $2^{\rm{nd}}$ \& $3^{\rm{rd}}$ derivatives and (2) uniqueness of its maximum in the region, behaves locally like a quadratic function around its maximum. In our case, $\mathcal{L}(\mathcal{S}, w)$ corresponds to the log-likelihood function $\ell\ell_{\mathcal{S}}(w)$ and so $w_0$ corresponds to a local maximum of $\ell\ell_{\mathcal{S}}(w)$. Since $\mathcal{L}(\mathcal{S}, w)$ is analytic and $w_0$ is the only local minimum of $\mathcal{L}(\mathcal{S}, w)$ in $\mathcal{M}(w_0)$, the training loss indeed can be considered locally quadratic.

Similar to \citet{langford2002not}, \citet{harvey2017nearly} and \citet{neyshabur2017exploring}, we apply the PAC-Bayes Theorem \citep{mcallester2003simplified} to derive a generalization bound for our metric. Specifically, we pick a uniform prior $\mathcal{P}$ over $w \in \mathcal{M}(w_0)$ according to the maximum entropy principle and pick the posterior $\mathcal{Q}$ of density $q(w) \propto e^{-|\mathcal{L}_0 - \mathcal{L}(\mathcal{S}, w)|}$ with $\mathcal{L}_0 \mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} \mathcal{L}(\mathcal{S}, w_0)$. Then Theorem \ref{main} bounds the expected generalization loss using $\gamma(w_0)$ (proved in Appendix B).

\begin{thm} \label{main}
Given $|\mathcal{S}| = N$, $\mathcal{D}$, $\mathcal{L}(\mathcal{S}, w)$ and $\mathcal{L}(\mathcal{D}, w)$ described in Sec. \ref{pform}, a local minimum $w_0$, the volume $V$ of $\mathcal{M}(w_0)$ sufficiently small, the Assumption \ref{assume_unique} \& \ref{assume} satisfied, and $\mathcal{P}, \mathcal{Q}$ defined above, for any $\delta \in (0, 1]$, we have with probability at least $1 - \delta$ that:
\begin{align*}
\mathop{\mathbb{E}}_{w \sim \mathcal{Q}} [\mathcal{L}(\mathcal{D}, w)] &\leq \mathop{\mathbb{E}}_{w \sim \mathcal{Q}} [\mathcal{L}(\mathcal{S}, w)] + 2 \sqrt{ \frac{2 \mathcal{L}_0 + 2 \mathcal{A} + \ln \frac{2N}{\delta}}{N - 1}} \\
\textrm{where}\ \mathcal{A} &= \frac{1}{4\pi e}W V^{\frac{2}{W}} \pi^{\frac{1}{W}} \exp\{\frac{\red{\gamma(w_0)}}{W}\}
\end{align*}

\end{thm}

Where $W$ is the number of model parameters (defined in Sec. \ref{pform}) and
$V$ the volume controlling the size of the neighborhood taken into account around $w_0$ (defined in Sec. \ref{metric}).
In short, Theorem \ref{main} shows that a lower $\gamma(w_0)$ indicates a local minimum $w_0$ of better generalization.

\subsection{Approximation} \label{appro}
As stated in Sec. \ref{fisher}, in practice an approximation of $\gamma(w_0)$ as $\widehat{\gamma}(w_0)$ is necessary, as calculating $\gamma(w_0)$ involves computing the product of all $W$ non-zero eigenvalues of the Fisher information matrix. Assume an imagined training set $\mathcal{S}'$ of size $W$ and a local minimum $w_0$ of $\mathcal{L}(\mathcal{S}', w)$; then $\ln |\mathcal{I}_{\mathcal{S'}}(w_0)|$ is well defined on the full-rank Fisher information denoted as $\mathcal{I}_{\mathcal{S'}}(w_0)$. In reality, we only have a training set $\mathcal{S} \subset \mathcal{S}'$ with $|\mathcal{S}|$ non-zero eigenvalues of the singular matrix $\mathcal{I}_{\mathcal{S}}(w_0)$.
Similar to the approach in \citet{karakida2019universal}, we propose to approximate eigenvalues of $\mathcal{I}_{\mathcal{S'}}(w_0)$ by the non-zero eignevalues of $\mathcal{I}_{\mathcal{S}}(w_0)$, or equivalently, as shown later, by the eigenvalues of sub-matrices of $\mathcal{I}_{\mathcal{S'}}(w_0)$.

First of all, we replace $y$ by its one-hot version $\tilde{y}$ defined in Sec. \ref{nota}, drastically reducing the cost of gradient calculation. This is reasonable since $y$ and $\tilde{y}$ are very close. With $\tilde{\mathcal{L}}(\mathcal{S}, w) \in \mathbb{R}^N$ defined in Sec. \ref{nota}, according to Equation \ref{eq:1}, we have $\mathcal{I}_{\mathcal{S'}}(w_0) \in \mathbb{R}^{W \times W}$ as:
\begin{align} \label{eq:one_hot_approx}
\mathcal{I}_{\mathcal{S'}}(w_0) &\approx \frac{1}{W}\sum_{(x, y) \in \mathcal{S}'} \nabla_{w}[\boldsymbol{\ell}_x(w_0)]_{\textbf{y}} \nabla_{w}[\boldsymbol{\ell}_x(w_0)]_{\textbf{y}}^T \nonumber \\
&\quad\ \textrm{where}\ \textbf{y} = \myargmax(y) \nonumber \\
&= \frac{1}{W} \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S'}, w)]^T \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S'}, w)] \nonumber \\
&= \frac{1}{W} \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S'}, w)]\ \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S'}, w)]^T
\end{align}
Let $\{\lambda_m\}_{m=1}^W$ denote the eigenvalues of $\mathcal{I}_{\mathcal{S'}}(w_0)$; then $\gamma(w_0) = \ln \prod_{m = 1}^W \lambda_m = \sum_{m = 1}^W \ln \lambda_m$. Ideally we want to perform a Monte-Carlo estimation of $\gamma(w_0)$ by randomly sampling $N' < N \ll W$ eigenvalues from $\{\lambda_m\}_{m=1}^W$, where $N$ is the size of $\mathcal{S}$. We denote the samples as $\{\lambda_n\}_{n=1}^{N'}$ and we have $\frac{W}{N'}\sum_{n=1}^{N'} \ln \lambda_n \approx \sum_{m=1}^W \ln \lambda_m$. Suppose the estimation is run $T$ times, we have $\lim_{T \rightarrow \infty} \frac{1}{T} \sum_{t=1}^T \frac{W}{N'}\sum_{n=1}^{N'} \ln \lambda_n = \gamma(w_0)$.

Then the eigenvalue approximation comes in. We sample $\mathcal{S}^t \subset \mathcal{S}$ i.i.d. with $|\mathcal{S}^t| = N'$ for $T$ times and define
\vspace{1pt}
\begin{equation}
\xi^t(w_0) \mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S}^t, w_0)] \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{S}^t, w_0)]^T \in \mathbb{R}^{N' \times N'}
\end{equation}
Notice that $\xi^t(w_0)$ is a principal sub-matrix of $W \mathcal{I}_{\mathcal{S'}}(w_0)$ by removing rows \& columns for data in $\mathcal{S}$ \textbackslash\ $\mathcal{S}^t$. According to Theorem 2 in Appendix C and properties of the spectral density of Fisher information \citep{pennington2018spectrum, sagun2017empirical, karakida2019universal}, one can well approximate the eigenvalues of $\mathcal{I}_{\mathcal{S'}}(w_0)$ by those of its sub-matrices. Therefore we
define the estimation $\widehat{\gamma}(w_0)$ as:
\begin{equation}
    \widehat{\gamma}(w_0) \mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} \frac{1}{T} \sum_{t = 1}^{T} \ln \big\vert\xi^t(w_0)\big\vert
\end{equation}
The relation between $\widehat{\gamma}(w_0)$ and $\gamma(w_0)$ is given as:
\begin{equation} \label{eq:2}
\quad \gamma(w_0) \approx \frac{W}{N'} \widehat{\gamma}(w_0) + W \ln \frac{1}{W}\ \ \textrm{as}\ \ T \rightarrow \infty \nonumber
\end{equation}

We leave the derivation of Equation \ref{eq:2} to Appendix C. In proposing $\widehat{\gamma}(w_0)$, we ignore the constants and irrelevant scaling factors. Empirically we find that given relatively large number of sample trials $T$, our metric $\widehat{\gamma}(\cdot)$ can effectively capture the generalizability of a local minimum even for a small $N'$ (details in Sec. \ref{exp_char} and Appendix D).

\section{Local Minima Regularization} \label{regularization}

Besides pragmatism, devising a practical regularizer based on $\gamma(w_0)$ also ``verifies'' our theoretical understanding of DNN training, helping the future improvement of the learning algorithms. Following the approximation scheme in Sec. \ref{appro}, it is natural to regularize $\gamma(w_0)$ during mini-batch learning by minimizing the product of $|\mathcal{B}|$ non-zero eigenvalues of the Fisher information computed $\mathcal{I}_{\mathcal{B}}(w_0)$, computed via the current batch $\mathcal{B}$, other than directly minimizing $\gamma(w_0)$.
However, this is far from practical due to the computation burden of:
\begin{enumerate}
    \item computing the eigenvalues in each training step
    \item computing second-order derivatives (i.e., computing the gradients of $\widehat{\gamma}(w_0)$ with respect to $w_0$)
\end{enumerate}
There is another major challenge. All of our theoretical analysis of $\gamma(\cdot)$ works on the grounds that the Assumption \ref{assume_unique} \& \ref{assume} are reasonable and satisfied, i.e., the largest $|\mathcal{B}|$ eigenvalues of $\mathcal{I}_{\mathcal{B}}(w_0)$ evaluated at the local minimum $w_0$ are non-zero. However, directly minimizing the product of these positive eigenvalues pays too much attention to the smallest eigenvalues, which can easily result in zero eigenvalues, raising singularity and thus violating the assumptions. Instead, we need the effort more spread out. A good choice is to minimize the trace norm $\textrm{tr}\big(\mathcal{I}_{\mathcal{B}}(w_0)\big)$, which provides an upper bound of the product of eigenvalues in the form of:
\[\prod_i \lambda_i \big(\mathcal{I}_{\mathcal{B}}(w_0)\big)^{1/|\mathcal{B}|} \leq \frac{1}{|\mathcal{B}|}\textrm{tr}\big(\mathcal{I}_{\mathcal{B}}(w_0)\big)\]
Although this bound is not be tight in general, we are tightening it when we minimize the trace norm. According to Corollary 1 in \citet{rodin2017variance}, we have:
\[
\frac{1}{|\mathcal{B}|}\textrm{tr}\big(\mathcal{I}_{\mathcal{B}}(w_0)\big) - \prod_i \lambda_i \big(\mathcal{I}_{\mathcal{B}}(w_0)\big)^{1/|\mathcal{B}|} \leq \sqrt{|\mathcal{B}| - 1}\ \sigma
\]
Where $\sigma$ denotes the standard deviation of the eigenvalues of $\mathcal{I}_{\mathcal{B}}(w_0)$. As pointed out in \citet{pennington2018spectrum, sagun2017empirical, karakida2019universal}, these eigenvalues are highly concentrated with only a few very large ``outliers'' which contribute the most to the variance. When we minimize the trace norm, i.e. the L1 norm of the eigenvalues, the largest few eigenvalues bear the most weight before they are reduced to a level that has the bound effectively tightened. Furthermore, computing the trace norm does not require computing eigenvalues; thus optimizing them removes the first computation burden.

Similar to the approach in Equation \ref{eq:one_hot_approx}, we approximate $y$ by its one-hot version $\tilde{y}$.
For simplicity, for the rest of this section, we denote $\tilde{y}$ as $y$ and correspondingly denote $\mathbb{E}_j [\tilde{\mathcal{L}}(\cdot, w_0)]_j$ as $\mathcal{L}(\cdot, w_0)$, where $\tilde{\mathcal{L}}$ is the one-hot vectorized loss defined in Sec. \ref{nota}.

Given a vector $x$, we have $\textrm{tr}(xx^T) = ||x||^2_2$. Therefore, we choose to approximate the trace norm as:
\begin{align*}
\textrm{tr}\big(\mathcal{I}_{\mathcal{B}}(w_0)\big) &\approx \frac{1}{|\mathcal{B}|} \sum_{(x, y) \in \mathcal{B}} \big\Vert\nabla_w \ell(f_{w_0}(x), y)\big\Vert^2_2
\end{align*}
To compute such quantity we need gradients for each individual data point.
We simplify this computation by grouping data points into batches and computing averaged gradients instead.
We randomly split $\mathcal{B}$ into $M$ sub-batches of equal size, namely $\{\mathcal{B}_i\}_{i=1}^{M}$. We define
\[\boldsymbol{g}_i \mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} \nabla_w \mathcal{L}(\mathcal{B}_i, w_0)\]
and then choose to optimize $\frac{1}{M}\sum_i^M \Vert\boldsymbol{g}_i\Vert^2_2$ instead of optimizing $\frac{1}{|\mathcal{B}|}\sum_{(x,y)\in\mathcal{B}}\Vert \nabla_w \ell(f_{w_0}(x), y) \Vert^2_2$, which drastically boosts the speed performance.

We deal with the second computation burden by adopting first order approximation. For any $w$, with a sufficiently small $\alpha > 0$, we have:
\[\tilde{\mathcal{L}}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i) \approx \tilde{\mathcal{L}}(\mathcal{B}_i, w) - \mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{B}_i, w)]\ \alpha \boldsymbol{g}_i\]
Thereby, we can estimate $\big\Vert\boldsymbol{g}_i\big\Vert^2_2$ by:
\begin{align*}
\alpha \big\Vert\boldsymbol{g}_i\big\Vert^2_2 &=  \frac{1}{|\mathcal{B}_i|}\sum_{j=1}^{|\mathcal{B}_i|}\Big[\mathbf{J}_{w}[\tilde{\mathcal{L}}(\mathcal{B}_i, w)]\ \alpha \boldsymbol{g}_i\Big]_j  \\
&\approx \frac{1}{|\mathcal{B}_i|}\sum_{j=1}^{|\mathcal{B}_i|}[\tilde{\mathcal{L}}(\mathcal{B}_i, w) - \tilde{\mathcal{L}}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i)]_j \\
&= \mathcal{L}(\mathcal{B}_i, w) - \mathcal{L}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i) \\
\end{align*}
Therefore, we propose to optimize the following regularized training objective for each mini-batch gradient descent step:
\begin{align} \label{eq:reg}
&\quad \mathcal{L}(\mathcal{B}, w) + \beta \mathcal{R}_{\alpha}(w)\ \ \textrm{where} \\ \mathcal{R}_{\alpha}(w) &\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}} \frac{1}{M} \sum_{i=1}^M \big[\mathcal{L}(\mathcal{B}_i, w) - \mathcal{L}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i)\big] \nonumber \\
&= \mathcal{L}(\mathcal{B}, w) - \frac{1}{M} \sum_{i=1}^M \mathcal{L}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i) \nonumber
\end{align}
Illustrated in Fig. \ref{fig:reg_illustrate}, an intuition is that Eq. \ref{eq:reg} penalizes a divergent set of gradients across samples in a mini-batch.

We omit any second order term when computing $\nabla_w \mathcal{R}_{\alpha}(w)$, simply by not back-propagating the gradient through $\boldsymbol{g}_i$.
We outline our regularized training step as Algorithm \ref{reg_algo}, which has 3 hyper-parameters: $\alpha$, $\beta$ and $M$.

\begin{algorithm}
\caption{Regularized Gradient Descent \footnotemark[1]} \label{reg_algo}
\begin{algorithmic}[1]
\Procedure{Update}{$w, \mathcal{B}$; $\alpha, \beta, M$}
\State $\{\mathcal{B}_i\}_{i=1}^M \gets \mathcal{B}$ \Comment{Split the mini-batch $\mathcal{B}$}
\For{$i \gets 1$ to $M$}
    \State $\boldsymbol{g}_i \gets  \nabla_w \mathcal{L}(\mathcal{B}_i, w_0)$
    \State $\boldsymbol{g}_i \gets \textrm{copy}(\boldsymbol{g}_i)$ \Comment{Stop the gradient\footnotemark[2]}
\EndFor
\State $\mathcal{R}_{\alpha}(w) \gets \frac{1}{M} \sum_{i=1}^M \big[\mathcal{L}(\mathcal{B}_i, w) - \mathcal{L}(\mathcal{B}_i, w - \alpha \boldsymbol{g}_i)\big]$
\State $\nabla_w \mathcal{L}_{\textrm{reg}} \gets \nabla_w [\mathcal{L}(\mathcal{B}, w) + \beta \mathcal{R}_{\alpha}(w)]$
\State Update weights $w$ with $\nabla_w \mathcal{L}_{\textrm{reg}}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\footnotetext[1]{Compatible with any gradient descent-based optimizer.}
\footnotetext[2]{Implemented as \verb|stop_gradient| in TensorFlow.}

\begin{figure}[ht]
  \vspace{-10px}
  \centering
  \includegraphics[width=\columnwidth]{reg.png}
   \vspace{-20px}
  \caption{An illustration of Algorithm \ref{reg_algo}. In essence, the regularizer guides the optimization process to areas with less divergent gradients of different data points within a mini-batch.}
  \label{fig:reg_illustrate}
\end{figure}

\section{Experiments} \label{exp}
We perform two sets of experiments to illustrate the effectiveness of our metric $\gamma(w_0)$. We demonstrate that: (1) the approximation $\widehat{\gamma}(w_0)$ captures the generalizability well across local minima; (2) our regularization technique based on $\gamma(w_0)$ provides consistent generalization gain for DNNs.

Throughout our theoretical analysis, we assume that label smoothing (LS) is applied during model training in order to obtain well-defined local minima (first mentioned in Sec. \ref{fisher}). In all our empirical evaluations, we perform both the version with LS applied and without. Results are very similar and so we stick to the version without LS to be consistent with the original setup in papers of the various DNNs that we used.
As a result, $\tilde{y}$ and $y$ refers to the same quantity.

\subsection{Experiments on Local Minima Characterization} \label{exp_char}

We perform comprehensive evaluations to compare our metric $\widehat{\gamma}(\cdot)$ with several others on ResNet-20 \citep{he2016deep} for the CIFAR-10 dataset (architecture details in Appendix E). Our metric consistently outperforms others in indicating local minima's generalizability. Specifically, \citet{sokolic2017robust} proposed a robustness-based metric used as a regularizer; \citet{wu2017towards} proposed to use Frobenius norm of the Hessian as a metric; \citet{keskar2016large} proposed a metric closely related to the spectral radius of Hessian. In summary, we compare 4 metrics, all evaluated at a local minimum $w$ given training set $\mathcal{S}$. All four metrics go for ``smaller values indicate better generalization''.
\begin{itemize}
  \item Robustness: $\frac{1}{N} \sum_{(x,y) \in \mathcal{S}} \big\Vert\textbf{J}_{x}[f_w(x)]\big\Vert^2_2$
  \item Frobenius norm: $\big\Vert \nabla^2_w \mathcal{L}(\mathcal{S}, w)\big\Vert_F^2$
  \item Spectral radius: $\rho( \nabla^2_w \mathcal{L}(\mathcal{S}, w))$
  \item Ours: $\widehat{\gamma}(w) =\frac{1}{T} \sum_{t = 1}^{T} \ln |\xi(\mathcal{S}^t, w_0)|$, $\ \mathcal{S}^t \subset \mathcal{S}$
\end{itemize}
Both the Frobenius norm and the spectral radius based metric are related to ours, as from Equation \ref{eq:b} we have $\big\Vert \nabla^2_w \mathcal{L}(\mathcal{S}, w)\big\Vert_F^2 = \big\Vert\mathcal{I}_{\mathcal{S}}(w)\big\Vert_F^2$ and $\rho(\nabla^2_w \mathcal{L}(\mathcal{S}, w)) = \rho(\mathcal{I}_{\mathcal{S}}(w))$. These two metric, however, are too expensive to compute for the entire training set $\mathcal{S}$; we instead calculate them by averaging the results for $T$ sampled $\mathcal{S}^t \subset \mathcal{S}$, similar to when we compute $\widehat{\gamma}(w)$. We leave details of how we exactly compute these metrics to Appendix D.

We perform evaluations in three scenarios, similar to \citet{neyshabur2017exploring, keskar2016large}. We compute the 4 metrics on different local minima arising due to (1) a confusion set of varying size in training, (2) different data augmentation schemes, and (3) different batch size.
\begin{itemize}
    \item In Scenario I, we randomly select a subset of 10000 images from CIFAR-10 as the training set and train the DNN with a confusion set consisting of images with random labels. We vary the size of the confusion set so that the resulting local minima generalize differently to the test set while all remain close-to-zero training losses. We consider confusion size of $0$, 1k, 2k, 3k, 4k and 5k. We calculate all metrics based on the sampled 10000 training images.
    \item In Scenario II, we vary the level of data augmentation. We apply horizontal flipping, denoted \code{flip-only}, random cropping from images with 1 pixel padded each side plus flipping, denoted \code{1-crop-f}, random cropping with 4 pixels padded each side plus flipping, denoted \code{4-crop-f} and no data augmentation at all, denoted \code{no-aug}. Under all schemes, the network achieves perfect training accuracy. All the metrics are computed on the un-augmented training set.
    \item In Scenario III, we vary the batch size. \citet{hoffer2017train} suggests that large batch sizes lead to poor generalization. We consider the batch sizes to be $128$, $256$, $512$ and $1024$.
\end{itemize}
The default values for the 3 variables are confusion size 0, \code{4-crop-f} and batch size 128. For each configuration in each scenario, we train 5 models and report results (average \& standard deviations) of all metrics as well as the test errors (in percentage). For the confusion set experiments, we sample a new training set and a new confusion set every time. In all scenarios, we train the model for 200 epochs with an initial learning rate 0.1, divided by 10 whenever the training loss plateaus. Within each scenario, we find the final training loss very small and very similar across different models and the training accuracy essentially equal to 1, indicating the convergence to local minima.

The results are in Figure \ref{fig:confusion}, \ref{fig:data_aug} and \ref{fig:bs} for Scenario I, II and III, respectively. Our metric significantly outperforms others and is very effective in capturing the generalization properties, i.e., a lower value of our metric consistently indicates a better generalizable local minimum.

\begin{figure*}[h!]
  \centering
  \includegraphics[width=17cm]{f1.PNG}
  \vspace{-5pt}
  \caption{Scenario I: Varied size of the confusion set. 5 models are trained for each size of the confusion set (x-axis). Solid lines are the average result; shaded areas represent the $\pm$ 1 standard deviation (same for Figure \ref{fig:data_aug} and \ref{fig:bs}). A larger confusion set leads to a higher test error, a trend well captured by our metric and the other two; the robustness based metric fails.}
  \label{fig:confusion}
\end{figure*}
\begin{figure*}[h!]
  \centering
  \includegraphics[width=17cm]{f3.PNG}
  \vspace{-5pt}
  \caption{Scenario II: Varied data augmentation schemes. Four different schemes are used. Our metric works well as an indicator of the test error while all the other metrics completely fail.}
  \label{fig:data_aug}
\end{figure*}
\begin{figure*}[h!]
  \centering
  \includegraphics[width=17cm]{f2.PNG}
  \vspace{-5pt}
  \caption{Scenario III: Larger batch size leads to worse generalization, captured by all the metrics except for the robustness based one.}
  \label{fig:bs}
\end{figure*}

\begin{table*}[h!]
  \caption{Test error (\%) on CIFAR-10/100. In general, a model with more parameters admits more space for regularization. The representation power of ResNet-20 is too limited for CIFAR-100 (resulting in poor convergence); so we ignore it in our experiments.}
  \vspace{7pt}
  \scriptsize
  \centering
  \setlength{\tabcolsep}{7pt}
  \begin{tabular}{lll|ll|ll|ll}
    \toprule
     & CNN & CNN+reg & WRN-28-2 & WRN-28-2+reg & DenseNet-k12 & DenseNet-k12+reg & ResNet-20 & ResNet-20+reg\\
    \midrule
     CIFAR-10 & 8.52 $\pm$ 0.23 & \textbf{7.55 $\pm$ 0.06} & 5.63 $\pm$ 0.20 & \textbf{5.15 $\pm$ 0.09} & 4.61 $\pm$ 0.08 & \textbf{4.37 $\pm$ 0.06} & 8.50 $\pm$ 0.31 & \textbf{7.89 $\pm$ 0.13}\\
     CIFAR-100 & 31.12 $\pm$ 0.35 & \textbf{29.27 $\pm$ 0.17} & 25.71 $\pm$ 0.24 & \textbf{23.88 $\pm$ 0.13} & 22.54 $\pm$ 0.32 & \textbf{22.23 $\pm$ 0.21} &  - & - \\
    \bottomrule
  \end{tabular}
 \label{tab:reg}
 \end{table*}

\begin{table*}[h!]
  \caption{Validation set error (\%) on $128\times 128$ down-sampled ImageNet classification. The better results are bolded.}
  \vspace{7pt}
  \scriptsize
  \centering
  \setlength{\tabcolsep}{7pt}
  \begin{tabular}{llll|llll}
    \toprule
    Top1 Error (\%) & Test & Train & Average Gap & Top5 Error (\%) & Test & Train & Average Gap \\
    \midrule
    WRN-18 & 35.52 $\pm$ 0.11 & 23.67 $\pm$ 2.05 & 11.85 & & 14.27 $\pm$ 0.02 & 7.33 $\pm$ 2.25 & 6.94 \\
    WRN-18+reg & \textbf{34.99 $\pm$ 0.10} & \textbf{24.0 $\pm$ 3.11} & \textbf{10.99} & & \textbf{13.85 $\pm$ 0.05} & \textbf{7.31 $\pm$ 1.07} & \textbf{6.54} \\
    \bottomrule
  \end{tabular}
 \label{tab:reg2}
 \end{table*}

\subsection{Experiments on Local Minima Regularization} \label{exp_reg}
We evaluate our regularizer on CIFAR-10, CIFAR-100 and the ImageNet classification task \citep{imagenet_cvpr09}. For CIFAR-10 \& CIFAR-100, we evaluate on four different network architectures including a plain CNN, ResNet-20, Wide ResNet \citep{wideresnet} and DenseNet \citep{huang2017densely}. We use WRN-28-2-B(3,3) from \citet{wideresnet} and the DenseNet-BC-k=12 from \citet{huang2017densely}.
We evaluate ImageNet classification on WRN-18-1.5 from \citet{wideresnet}.
In specific, we follow \citet{sokolic2017robust} to down-sample all images to $128\times 128$ and apply standard data augmentations.
See Appendix E for architecture and training details.
We denote the four networks as CNN, ResNet-20, WRN-28-2 / WRN-18 and DenseNet-k12, respectively.

For the three hyper-parameters $\alpha, \beta, M$ in our proposed Algorithm \ref{reg_algo}, we find $\alpha$ and $M$ quite robust and manually set $\alpha = 0.0001$, $M = 8$ in all experiments and select $\beta$ by validation via a 45k/5k training data split for each of the network architecture \& dataset pair. In specific, we consider $\beta \in \{1, 5, 10, 20, 30, 40, 50, 75, 100\}$. We keep all the other training hyper-parameters, schemes as well as the setup identical to those in their original paper whenever possible (details in Appendix E). We train 5 separate models for each network-dataset combination on CIFAR-10 and CIFAR-100 and train 3 models for ImageNet. We report the test errors in percentage (mean $\pm$ std.) in Table \ref{tab:reg} and \ref{tab:reg2}, where ``+reg'' indicates training with our regularizer applied. The results demonstrate that our method provides consistent generalization improvement for a wide range of DNNs.

\subsubsection{Time Complexity for Algorithm \ref{reg_algo}}
We benchmark WRN-18 on the down-sampled ImageNet classification dataset with 2 Nvidia 2080 Ti GPUs and a batch size of 128. With parallelization, the average training time per mini-batch is 185.7ms without regularizer applied vs. 285.6ms with regularizer applied.
It only takes around 1.5x longer time per gradient update for Algorithm \ref{reg_algo}.

By ablation study, we find that our regularizer works the best in the mid and late stage of DNN training, e.g., we only use the regularized update after the first learning rate drop in all of our experiments. In the beginning stage where the optimization process is not stable, our regularizer can result in great numerical errors. By only applying Algorithm \ref{reg_algo} during the later stages, the training speed can be further increased by a large margin.

\subsubsection{The Choice of The Optimizer}
As described in Algorithm \ref{reg_algo}, our proposed regularizer is not tied to a specific optimizer. We perform experiments with SGD+Momentum because it is chosen to be used in ResNet, WRN, and DenseNet, helping all of them achieve current or previous state-of-the-art results. Our regularizer aims to find better ``flatter'' minima to improve generalization whereas adaptive optimization methods such as Adam \citep{kingma2014adam} and AdaGrad \citep{duchi2011adaptive} try to boost up convergence, yet usually at the cost of generalizability. Recent works \citep{wilson2017marginal, keskar2017improving} show that adaptive methods generalize worse than SGD+Momentum. In specific, very similar to our setup, \citet{keskar2017improving} demonstrates that SGD+Momentum consistently outperforms the others on ResNet and DenseNet for CIFAR-10 and CIFAR-100. Other approaches that also utilize local curvature to improve SGD, such as the Entropy-SGD \citep{chaudhari2016entropy} mentioned in Sec. \ref{related_work}, have empirical results rather preliminary compared to ours.

\begin{table}[h!]
  \caption{The proposed metric computed on local minima obtained with or without applying the proposed regularizer. Each entry represents mean $\pm$ std. among 5 runs. Smaller values are bolded.}
  \vspace{7pt}
  \scriptsize
  \centering
  \setlength{\tabcolsep}{4pt}
  \begin{tabular}{l|l|l|l}
    \toprule
     & ResNet-20 & WRN-28-2 & DenseNet-k12 \\
    \midrule
     w/o reg. & -979.3 $\pm$ 22.3 & -689.6 $\pm$ 24.9 & -850.3 $\pm$ 23.5 \\
     with reg. & \textbf{-1138.1 $\pm$ 11.0} & \textbf{-748.7 $\pm$ 21.3} & \textbf{-886.2 $\pm$ 20.5} \\
    \bottomrule
  \end{tabular}
 \label{tab:metric}
 \end{table}

\subsubsection{Generalization Boost As a Result of Better Local Minima}
We perform a sanity check to illustrate that our regularizer indeed induces better local minima characterized by our metric, i.e., our proposed regularizer is consistent with our proposed metric. For ResNet, Wide-ResNet and DenseNet trained on CIFAR-10, we compute the metric on local minima obtained with or without applying the regularizer.
In specific, our regularizer has an impact on the optimization process, leaving training loss slightly different for models with or without the regularizer.
To ensure our assumption that those local minima have similar close-to-zero training loss, before computing $\hat{\gamma}$ for each model, we normalize and scale the softmax output for each individual training sample. This operation makes comparison between different DNN models robust without changing their underlying behaviors.
Table \ref{tab:metric} shows that the resulting generalization boost aligns with what captured by our metric.

\section{Conclusion and Future Work}
In this paper, we show a bridge between the field of deep learning theory and regularization methods with respect to the generalizability of local minima. We propose a metric that captures the generalization properties of different local minima and provide its theoretical analysis including a generalization bound. We further derive an efficient approximation of the metric and a practical and effective regularizer. Empirical results demonstrate our success in both capturing and improving the generalizability of DNNs.

Moreover, we find that our proposed regularizer might be further simplified and a dynamic scheduling of the hyper-parameter $\beta$ can provide even more improvement to the generalization performance.
In general, our exploration promises a direction for future work on the regularization and optimization of DNNs.

\paragraph{Acknowledgment}
This work was supported in part by NSF awards CNS-1730158.


\section*{}
\vspace{-1cm}

\footnotetext{\textit{Nanoscience Center, Department of Physics, University of Jyv\"askyl\"a, 40014 Jyv\"askyl\"a, Finland. E-mail: pekka.koskinen@iki.fi}}

\footnotetext{\dag~Electronic supplementary information (ESI) available. See DOI: 10.1039/c9nr08533e

Mainstream two-dimensional (2D) materials research has long focused on layered van der Waals materials~\cite{miro14,lin16,novoselov19}. These materials have rigid covalent in-plane bonding and weak van der Waals out-of-plane bonding. They are interesting due to simple fabrication~\cite{novoselov05}, intrinsic properties~\cite{geim09,novoselov16}, possibility of materials design~\cite{geim13,Cao2018}, and a myriad of applications~\cite{deng16,fiori14,xu13}. However, yet a nascent practice is to use these 2D materials as templates to stabilize other, more exotic nanostructures.

Stabilization occurs particularly well in pores, which are ubiquitous in 2D materials~\cite{russo12}. Pores enable stabilizing 2D nanostructures that otherwise would be unstable. Notably, pores in graphene have been used to stabilize even free-standing patches of 2D metals (Fig. \ref{fig:fig1}a)\cite{Wang2012a,zhao14,Ling2015a}. Here 2D metal patches mean monolayer clusters of elemental metals that are characterized by metallic bonding, homogeneous electron density, and superatomic states~\cite{nevalaita18,lin_PRL_09,lin_PRB_10,Stiehler2013}. As elemental metals usually prefer 3D clustering~\cite{koskinen07}, free-standing 2D metal patches represent materials under conditions so extraordinary that they are bound for investigations, as their properties and applications could differ markedly from the ones found in supported metal monolayers~\cite{yin09,ling15,chen18, ma18}.

However, synthesis and experimental control of 2D metal patches is challenging and still immature. For example, the usual mechanical exfoliation of 2D metals is unviable~\cite{coleman11,nicolosi13,Tang2014a}. Only few experiments have demonstrated inklings of success, reported as small graphene-stabilized patches of Fe~\cite{zhao14}, MoS$_2$-stabilized patches of Mo~\cite{zhao18}, as well as somewhat related graphene-stabilized patches of and Zn and Cu oxides~\cite{quang15, yin17}. Consequently, present research is best driven forward by simulations and modeling. Nevertheless, so far simulations have been limited to scattered elements and to effectively infinite membranes~\cite{yang15, yang15b, yang16, koskinen15, antikainen17,nevalaita18,nevalaita18b}. As a result, a practical and coherent understanding of the energetic trends in the patches as well as the most promising candidates among elemental metals remain unknown.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{fig1.pdf}
    \caption{Elemental 2D metal patches in covalent templates. a) Schematic of an atomically thin, free-standing layer of close-packed elemental 2D metal monolayer patching a pore in graphene. b) The 45 elemental metal candidates explored in this work.}
    \label{fig:fig1}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{fig2.pdf}
\caption{Relative stability of 2D and 3D clusters in gas-phase, without the presence of stabilizing pores.  a) Energy per atom for 2D and 3D clusters of Cu according to the liquid drop model [eqs~(\ref{eq:E3D}) and (\ref{eq:E2D})]. b) Same as panel a for Ir. c) Dimensionless parameter $\Delta$ [eqn~(\ref{eq:delta})] visualized across the periodic table; smaller value stands for greater intrinsic 2D stability.}
    \label{fig:fig2}
\end{figure}

Therefore, in this work we address the fundamental question: \emph{Which elemental metals enable the most stable 2D patches in the pores of covalent 2D materials?} Using density-functional theory (DFT), we investigate finite patches of 2D metals and their interaction with a graphene template. The trends in key properties are explored across the periodic table, among 45 elemental metals (Fig.~\ref{fig:fig1}b). By combining DFT calculations with a liquid drop model, we are able to identify parts of the periodic table that hold the most prominent 2D metals. It turns out that, for graphene pores, the largest and most stable 2D metal patches are provided by Zn, Ag, and Au, the best candidate being Cu.

To investigate the energetic stability trends of 2D metal patches, we resort to the reliable liquid drop model~\cite{perdew91}. We begin by analyzing the stabilities of 2D and 3D clusters first in gas-phase, without the template. According to the model, the formation energy of a 3D cluster with $N$ atoms is
\begin{equation}
E(N) = \alpha V(N)+\sigma A(N) + E_c(N),
\end{equation}
where $V$ is the volume, $\alpha$ is the bulk energy density, $A$ is the surface area, $\sigma$ is the surface energy density, and $E_c$ is energy due to surface curvature. For simplicity, we choose closely packed, spherical clusters with radius $R$ and neglect the curvature energy term. Thus, the formation energy per atom for an $N$-atom cluster becomes
\begin{equation}
\varepsilon_\text{3D}^\text{gas}(N) = - \varepsilon^\text{3D}_\text{coh}+\sigma A(N)/N,
\label{eq:E3D}
\end{equation}
where $\varepsilon^\text{3D}_\text{coh}$ is the 3D bulk cohesion, $A(N)=c_\text{3D}d_\text{3D}^2 N^{2/3}$, $c_\text{3D}=(18\pi)^{1/3}\approx 3.84$, and $d_\text{3D}$ is the 3D bond length. Analogously, liquid drop model for the formation energy of a disc-shaped 2D patch is
\begin{equation}
\varepsilon_\text{2D}^\text{gas}(N) = - \varepsilon^\text{2D}_\text{coh}+\lambda L(N)/N,
\label{eq:E2D}
\end{equation}
where $\varepsilon^\text{2D}_\text{coh}$ is the 2D bulk cohesion, $\lambda$ is the edge energy, $L(N)=c_\text{2D}d_\text{2D}N^{1/2}$ is the length of the edge, $c_\text{2D}=(12 \pi^2)^{1/4}\approx 3.30$, and $d_\text{2D}$ is the 2D bond length.

Since the 3D data for eqn~(\ref{eq:E3D}) can readily be found from the literature and the 2D data for eqn~(\ref{eq:E2D}) is found from our recent work \cite{nevalaita18,nevalaita18b}, we may proceed to investigate the relative stability of 2D and 3D gas-phase clusters quantitatively. Clearly, 3D clusters are more stable than 2D clusters and their relative stability increases faster upon increasing $N$ (Figs.~\ref{fig:fig2}a and \ref{fig:fig2}b). For some metals, such as Cu, crossover from 2D to 3D occurs at $N$ too small to have any quantitative meaning (Fig.~\ref{fig:fig2}a), while for most metals, such as Ir, crossover is absent altogether and the ground state is 3D for all $N$~(Fig.~\ref{fig:fig2}b). Not surprising, these trends are in accordance with the long known prevalence of 3D ground states in metal clusters~\cite{hakkinen02}.

However, despite the predictability of the 2D/3D energy ordering, it is illustrating to equate (\ref{eq:E3D}) and (\ref{eq:E2D}) to obtain an expression for the 2D/3D crossover. That expression
\begin{equation}
(\epsilon^\text{3D}_\text{coh}-\epsilon^\text{2D}_\text{coh})-\sigma c_\text{3D}d^2_\text{3D}N^{-1/3}+\lambda c_\text{2D} d_\text{2D}N^{-1/2}=0
\label{eq:cross}
\end{equation}
is simplified by two approximations. First, our previous work shows that for close-packed 2D metals $d_\text{3D}\approx d_\text{2D}$~\cite{nevalaita18}. Second, because close-packed structures are nearly isotropic, the surface energy and edge energy are related by $\sigma \approx \lambda/d_\text{2D}$. (See ESI\dag~for derivation of this relation.) Thus, the 2D/3D crossover occurs when
\begin{equation}
\Delta -c_\text{3D}N^{-1/3}+c_\text{2D}N^{-1/2}=0,
\label{eq:crossing}
\end{equation}
where
\begin{equation}
\Delta = (\varepsilon^\text{3D}_\text{coh}-\varepsilon^\text{2D}_\text{coh})/(d_\text{2D}\lambda)
\label{eq:delta}
\end{equation}
is a dimensionless parameter, the smallness of which characterizes the intrinsic relative stability of 2D patches.

The crossover occurs if Eq. (5) has a positive real solution, requiring $\Delta < 4/\sqrt{27}\approx 0.77$, which nominally holds only for few elements. Nevertheless, visualizing $\Delta$ across the periodic table gives us insight into elements' \emph{intrinsic} 2D stabilities (Fig. \ref{fig:fig2}c). As an intermediate result, most of the prominent candidates for elemental 2D metals are located at the end of the transition metal series. This trend is corroborated by the known disposition of Au towards 2D ground state geometries~\cite{hakkinen02,koskinen_NJP_06,koskinen07}.

The gas-phase results above provide a good starting reference, but our focus lies in the stabilization due to the interaction with pore edges. For concreteness, we will restrict ourselves to graphene pores~\cite{russo12}, because of graphene's mechanical rigidity and the feasibility of controlling pore sizes~\cite{Robertson2015,Wang2014a,kudin_PRB_01,koskinen_PRB_10b}.

The stabilization of 2D clusters occurs only when the interaction between the metal atoms and the pore edge is exothermic and suitably directional. This requirement is fulfilled by graphene, which prefers in-plane adsorption at edges~\cite{Malola2009,antikainen17} as opposed to on-top adsorption~\cite{pastewka_JPS_13}. On-top adsorption energies, the values of which depend largely on the DFT functional used~\cite{pasti18}, remain below $2$~eV, mostly even below $0.2$~eV, whereas edge adsorption energies can be up to $7$~eV, depending on edge type and adsorption site (Figs.~\ref{fig:fig3}a and Fig.~S3; see ESI\dag~for computational details). These trends indicate that the nature of metal-graphene interaction at pore edges is well-suited for stabilizing 2D patches.

A caveat, however, is that the interaction and in-plane adsorption must not be too strong. Too strong and inapt interaction may cause graphene to swallow metal atoms altogether and give rise to metal carbides. We must therefore extend our investigations to various metal carbides and exclude elements that prefer carbide formation. To do this, we define carbide formation energy
\begin{equation}
    \Delta \varepsilon_\text{carbide}=[E_\text{carbide}-(N_C\mu_\text{gr}+N_M\mu_M)]/N_M,
    \label{eq:ecarb}
\end{equation}
where $E_\text{carbide}$ is the total energy of a carbide with $N_M$ metal atoms and $N_C$ C atoms, $\mu_\text{gr}$ is the chemical potential of C in graphene, and $\mu_M$ is the chemical potential of metal atoms in 2D bulk. $\Delta \varepsilon_\text{carbide}$ was calculated for all metals with three different carbide geometries (Fig.~S5\dag). Any occurrence of negative carbide formation energy implies that 2D metal-graphene interface is unstable; this consideration eliminates many of our candidates (Fig.~\ref{fig:fig3}b).

Above we considered edge adsorption of single atoms, which is different from the interaction between graphene and a finite 2D metal patch. To prevent interface buckling and the concomitant structural out-of-plane perturbations, interface should be commensurate whereby lattice mismatch between graphene and the metal patch should be minimal. A look at the matching of lattice constants thus provides yet another view into the most promising candidates for stable 2D patches with graphene (Fig.~\ref{fig:fig3}c). Incidentally, the picture of promising elements looks nearly the same for hexagonal boron nitride template, which has only $1.8$~\%\ lattice mismatch with graphene~\cite{Woods2014}.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.9\columnwidth]{fig3.pdf}
    \caption{Viewpoints to interaction between elemental metals and graphene. a) Metal atom adsorption energies at graphene edges. The values shown are the minima among adsorption energies at hollow and bridge positions of armchair edges (Fig.~S3\dag). b) Carbide formation energies [eqn~(\ref{eq:ecarb})]. The values shown correspond to the minima of formation energies among three carbide configurations (Fig.~S5\dag). c) Absolute strains in 2D metals due to lattice mismatch relative to graphene. d) Interface energies between selected small-strain elemental 2D metals and graphene. The values shown are the minima among energies of four different metal-graphene interface configurations~(Fig.~S4\dag). }
    \label{fig:fig3}
\end{figure}

Consequently, we chose nine metals with best-matching lattice constants for further analysis. Although V had the smallest lattice mismatch, it was dropped because of its negative carbide formation energy (Fig.~\ref{fig:fig3}b). We modeled the interface between metal and graphene using a ribbon model with different interface configurations (Fig.~S2\dag). Metal-graphene interface energy was defined as
\begin{equation}
\lambda_\text{if} = (E_M+E_\text{gr}-E_\text{M+gr})/L_\text{if},
\end{equation}
where $E_M$ is the energy of the 2D metal ribbon, $E_\text{gr}$ is the energy of the graphene ribbon, $E_\text{M+gr}$ is the energy of fused ribbons, and $L_\text{if}$ is the length of the ribbon interface. Note that the interface energy entails the strain energy of the 2D metal. All of the resulting interface energies are positive, which indicates energetically stable metal-graphene interfaces (Fig.~\ref{fig:fig3}d); the dependence on the precise microscopic configuration was moderate (Fig.~S4\dag). Contrary to the intrinsic stability of gas-phase clusters (Fig.~\ref{fig:fig2}c), interface energies favor elements around the middle of the transition metal series.

\begin{figure}[ht]
    \centering
    \includegraphics[width=\columnwidth]{fig4.pdf}
    \caption{Relative stability of 2D and 3D clusters interacting with a graphene pore. a) Different interaction geometries: 3D clusters with on-top ($E_\text{ads}^{1}$) and pore edge ($E_\text{ads}^{2}$) adsorption, and 2D patches inside pores ($E_\text{ads}^{3}=\lambda_\text{if}L$).
    b) Energy per atom for the three geometries for Cu ($E_\text{ads}^{1}=1$~eV, $E_\text{ads}^{2}=15$~eV, and $\lambda_\text{if}=0.62$~eV/\AA).
    c) Energy per atom for the three geometries for Ir ($E_\text{ads}^{1}=1$~eV, $E_\text{ads}^{2}=15$~eV, and $\lambda_\text{if}=0.94$~eV/\AA). The blue/red shading denotes the 2D/3D crossover.}
    \label{fig:fig4}
\end{figure}

We now have enough data to extend the liquid drop model to include interaction with graphene. This extension simply means the addition of cluster adsorption energies into eqs~(\ref{eq:E3D}) and (\ref{eq:E2D}), which yields the energy per atom
\begin{equation}
\varepsilon_\text{3D}(N) = - \varepsilon^\text{3D}_\text{coh}+\sigma A(N)/N + E_\text{ads}/N
\label{eq:E3D2}
\end{equation}
for 3D clusters and
\begin{equation}
\varepsilon_\text{2D}(N) = - \varepsilon^\text{2D}_\text{coh}+(\lambda-\lambda_\text{if})L(N)/N
\label{eq:E2D2}
\end{equation}
for 2D clusters, assuming that the edge and interface lengths are equal. Note that we assume a pre-existing pore and therefore omit energy contributions arising from the pore itself and that because the contribution from $E_\text{ads}$ scales as $\propto N^{-1}$, its influence at large $N$ is minor. Because $\lambda<\lambda_\text{if}$ for all metals (Fig.~\ref{fig:fig3}d), eqn~(\ref{eq:E2D2}) demonstrates explicitly the manner in which the pore stabilizes the patch with any finite $L$. When plotting the relative stability of 2D and 3D clusters using eqs~(\ref{eq:E3D2}) and (\ref{eq:E2D2}), the stabilizing effect of the pores becomes obvious (Figs.~\ref{fig:fig4}b and \ref{fig:fig4}c). For instance, 2D patches becomes energetically favored over 3D clusters with $N < 125$ for Cu and with $N < 45$ for Ir. The liquid drop model thus suggests that pores are able to stabilize 2D patches of respectable sizes.

\begin{figure}[!t]
    \centering
    \includegraphics[width=\columnwidth]{fig5.pdf}
    \caption{Identifying the best elements for stable 2D metal patches. a) Summative figures of merit of 2D stability for the 45 metal candidates. Smaller number means better 2D stability; numbers are rounded to closest integers.~(ESI\dag). b) Energy per atom for $N$-atom Cu clusters in three geometries: 3D clusters with on-top adsorption (red triangles; top panel on right), 3D clusters with pore edge adsorption (red squares; middle panel on right), and 2D patches with pore adsorption (blue circles; bottom panel on right). The lines represent the liquid drop model expressions (\ref{eq:E3D2}) and (\ref{eq:E2D2}) and symbols DFT values; the liquid drop model parameters are the same as in Fig.~\ref{fig:fig4}b. c) Snapshot views of a $37$-atom Cu patch in a graphene pore. d) Top and side views of Cu and C atom trajectories during a $36.9$~ps molecular dynamics simulation at $T=300$~K, demonstrating a kinetically stable patch. The snapshots in panel c are taken from this trajectory.}
    \label{fig:fig5}
\end{figure}

To answer the the fundamental question about the best candidates for stable 2D metal patches, let us then summarize our findings. We have investigated stability of 2D patches from four viewpoints: \emph{i}) intrinsic relative stability of the 2D metal (Fig.~\ref{fig:fig2}c), \emph{ii}) commensurability and the resulting strain at the interface (Fig.~\ref{fig:fig3}c), \emph{iii}) possibility of carbide formation (Fig.~\ref{fig:fig3}b), and \emph{iv}) the chemical bonding at the interface (Fig.~\ref{fig:fig3}a and \ref{fig:fig3}d). Because some viewpoints go beyond pure energetics, a balanced ranking of the elements must consider them simultaneously. We therefore united all the viewpoints into a single summative figure of merit, and ranked the metals accordingly (see ESI\dag~for description). According to the figure of merit, the foremost metals for stable patches are Zn, Ag, and Au, but the best one is Cu (Fig.~\ref{fig:fig5}a).

Adopting Cu as the best candidate, let us first go backwards to make a necessary validation of the extended liquid drop model. We created $27$ Cu clusters with $N=5\ldots 145$ and calculated the formation energies explicitly by DFT (ESI\dag). While the clusters do not
necessarily represent global energy minima, the energy differences between low-energy isomers are small enough (< 15 meV/atom) to have any influence on the energetic trends~\cite{grigoryan06,assadollahzadeh08,zhang19}. As before (Fig.~\ref{fig:fig4}a), we considered three different geometries: $1$) 3D clusters with on-top adsorption, $2$) 3D clusters with edge adsorption, and $3$) 2D patches in graphene pores (ESI\dag). The extended liquid drop model seems to capture the trends in the energetics reasonably well, even upon assuming an $N$-independent 3D adsorption energy~(Fig.~\ref{fig:fig5}b). Most important, DFT predicts energetically stable 2D patches up to $N\approx 120$, corresponding to an area of $\sim 8$~nm$^2$, and supports the predictions of the liquid drop model. While there are variations in the DFT energies due to differences in atomic arrangements, the fair agreement in Fig.~\ref{fig:fig5}b suggests that the trends are correct and indicates that stability could be achievable for respectable patch sizes; besides, the stability would even improve upon making the patches oblong.

In addition to pure energetics, we made preliminary investigation of kinetic stability. We conducted molecular dynamics simulation of an $N=37$ Cu patch for $36.9$~ps using a Langevin thermostat at $T=300$~K. In addition to corroborating previous simulations that have demonstrated kinetically stable infinite 2D Cu membrane~\cite{yang16}, our simulations suggest kinetically stable graphene-metal interface as well (see Figs.~\ref{fig:fig5}c and \ref{fig:fig5}d as well as Supplemental Movie ESI\dag).

To conclude, we have used liquid drop model in conjunction with DFT to show that pores in covalent 2D templates such as graphene could be used to stabilize 2D metal patches in respectable sizes.
Although the best candidate, Cu, is familiar from graphene synthesis~\cite{Li2009b}, its prominence in this context is less clear. Anyhow, given the demonstrated usefulness of the liquid drop model, the results are swiftly generalizable to other covalent templates. The availability of free-standing, finite metal patches could serve as a platform for a number of applications including catalysis, sensing, bioimaging, photothermal therapy, solar cells, and electrical contacting~\cite{chen18}. Although we have predicted that the eventual patch structures ought to be stable, addressing the actual, experimental route to ultimately reach them is beyond our scope. However, pursuing this route resolutely is just what we legitimately may now propose.

\section*{Conflicts of interest}
There are no conflicts to declare.

\section*{Acknowledgements}
We acknowledge the Academy of Finland for funding (Project 297115).

\balance

\bibliographystyle{rsc}


\section{Introduction}

Let $a$ and $n$ be two relatively prime integers. The index of $a$ modulo $n$, denoted by $ind_{n}(a)$, is the smallest positive number $k$ such that $a^{k}\equiv 1\pmod{n}$. If $ind_{n}(a)=\phi(n)$, where $\phi$ is the Euler's totient function, then we say $a$ is a primitive root modulo $n$. Primitive roots are widely used in cryptography, such as attacking the discrete log problem, key exchange problem and many other public key cryptosystem. The primitive root theorem identifies all moduli which primitive roots exist, that is,
$$1,\ 2,\ 4,\ p^{\alpha}\ \text{and }2p^{\alpha}$$
where $p$ is an odd prime and $\alpha$ is a positive integer. For an odd prime $p>3$, Gauss obtained the following congruences in Article $81$,

\begin{theorem}[Gauss]\label{gausspr1}
$$\prod_{g\text{ is a primitive root mod }p}g\equiv1\pmod p.$$
\end{theorem}

\begin{theorem}[Gauss]\label{gausspr2}
$$\sum_{g\text{ is a primitive root mod }p}g\equiv\mu(p-1)\pmod p.$$
\end{theorem}

Few papers have been published concerning these two awesome congruences until Cai \cite{cai2017book} recently proved

\begin{theorem}[Cai]\label{caipr1}
Let $p$ be a prime greater than $3$ and $\alpha$ a positive integer. Then
$$\prod_{g\text{ is a primitive root mod }p^{\alpha}}g\equiv1\pmod{p^{\alpha}}.$$
\end{theorem}

\begin{theorem}[Cai]\label{caipr2}
Let $p$ be an odd prime and $\alpha$ a positive integer. Then
$$\sum_{g\text{ is a primitive root mod }p^{\alpha}}g\equiv\mu(p-1)\phi(p^{\alpha-1})\pmod{p^{\alpha}}.$$
\end{theorem}

We will continue Cai's work and give many interesting congruences concerning the sum and product of integers with same index. In this article, $\alpha$ and $\beta$ represent positive integers; $m$ is some given positive intger; $U_m$ is the set of invertible elements in $\mathbb{Z}/m\mathbb{Z}$; $\lambda(m)$ is the Carmichael function, namely, $\lambda(m)=$ min$\{k>0:$ $a^k\equiv1\pmod m$ for any $a\in U_m$ $\};$ $\delta$ is a positive divisor of $\lambda(m)$.

It is not a tough work to obtain the product version as follows.

\begin{theorem}\label{th:2.3.1}
\begin{equation}\label{eq:2.3.1}
\prod_{\substack{a\in U_m\\ind_m(a)=\delta}}a\equiv\begin{cases}-1\pmod m,&\text{ If }\delta=2\text{ and }m\text{ has a primitive root,}\\1\pmod m,&\text{ otherwise.}\end{cases}
\end{equation}
\end{theorem}

As for the sum, it is not a simple congruence written in one line. However, for fixed $\delta$, one can easily get the following results.

\begin{theorem}\label{th:2.3.2}
(1) $$\sum_{\substack{a\in U_m\\ind_m(a)=1}}a\equiv 1\pmod m;$$
(2) $$\sum_{\substack{a\in U_m\\ind_m(a)=2}}a\equiv -1\pmod m;$$
(3) $$\sum_{\substack{a\in U_m\\ind_m(a)=\delta\\4\mid \delta}}a\equiv 0\pmod m.$$
\end{theorem}

By applying this theorem, we obtain a result related to Fermat primes.

\begin{corollary}
If the regular $m$-gon can be constructed by compass and straightedge, then

$$\sum_{\substack{a\in U_{m}\\ind_{m}(a)=\delta}}a\equiv 0\pmod m.$$

\end{corollary}

To give the congruence mod general integer $m$, we first list the ones with some special moduli. For integers with primitive roots, we have

\begin{theorem}\label{th:2.3.3}
(1) $$\sum_{\substack{a\in U_m\\ind_2(a)=\delta}}a\equiv 1\pmod 2;$$
(2) $$\sum_{\substack{a\in U_m\\ind_4(a)=\delta}}a\equiv (-1)^{\delta+1}\pmod 4;$$
(3) $$\sum_{\substack{a\in U_m\\ind_8(a)=\delta}}a\equiv (-1)^{\delta+1}\pmod 8;$$
(4) If $\alpha>3$, $\delta=2^{\beta}$ and $1<\beta\leq\alpha-2$. Then
$$\sum_{\substack{a\in U_{2^{\alpha}}\\ind_{2^{\alpha}}(a)=\delta}}a\equiv 0\pmod{2^{\alpha}};$$
(5) Let $p$ be an odd prime. Then
$$\sum_{\substack{a\in U_{p^{\alpha}}\\ind_{p^{\alpha}}(a)=\delta}}a\equiv \mu((\delta,p-1))\phi(p^{ord_{p}(\delta)})\pmod{p^{\alpha}}.$$
where $ord_{p}(n)=max\{k:p^k\mid n\}.$
\end{theorem}

For a product of two distinct odd primes, we have

\begin{theorem}\label{th:2.3.4}
Let $p$ and $q$ be two distinct odd primes. Then
$$\sum_{\substack{a\in U_{pq}\\ind_{pq}(a)=\delta}}a\equiv \mu((\delta,p-1))I((\frac{\delta}{(\delta,p-1)},p-1))\frac{\phi(\delta)}{\phi((\delta,p-1))}\pmod p.$$
\end{theorem}

We end this section with the following congruence.

\begin{theorem}\label{th:2.3.5}

If $m=\prod_{i=1}^{k}p^{\alpha_i}$, then for any $p^{\alpha}\parallel m$,
\begin{align}\label{eq:th:2.3.5.2}
\sum_{\substack{a\in U_{m}\\ind_{m}(a)=\delta}}a
\equiv\mu((\delta,p-1))I((\frac{\delta}{(\delta,p-1)},p-1))F(\phi(p_1^{\alpha_1}),\cdots,\phi(p_k^{\alpha_k})&;\frac{\delta}{(\delta,p-1)})\\\nonumber&\pmod{p^{\alpha}}.
\end{align}
where $F(a_1,\cdots,a_k;n)=\prod_{p^{\alpha}\parallel n}\prod_{ord_p(a_i)<\alpha}p^{ord_{p}(a_i)}\prod_{ord_p(a_i)\geq\alpha}\phi(p^{\alpha})$.
\end{theorem}

Though this result seems complex, it shows that there is no need to compute the sum of integers with same index modulo $m$ by finding each of them.

\section{Preliminaries}

The following two lemmas are well known facts in the study of primitive roots. (See \cite{apostol}.)

\begin{lemma}\label{apostolind1}
Let $a$ belong to $U_{m}$. Then for any positive integer $k$,
$$ind_m(a^k)=\frac{ind_{m}(a)}{(k,ind_{m}(a))}.$$
\end{lemma}

\begin{lemma}\label{apostolind2}
Let $m$, $n$ and $a$ be positive integers, relatively prime in pairs. Then $ind_{mn}(a)=[ind_{m}(a),ind_{n}(a)]$. In particular, $\lambda(mn)=[\lambda(m),\lambda(n)],$ where $[m,n]$ represents the least common multiple of $n$ and $m$.
\end{lemma}

The Dirichlet convolution defined on the set of arithmetical functions: $\mathcal{A}:=\{f:\mathbb{Z}^{+}\to\mathbb{C}\}$ is very important in number theory, especially in analytic number theory. Two arithmetical functions $f$ and $g$ are combined by the convolution as follows
$$(f*g)(n)=\sum_{ab=n}f(a)g(b).$$
All the functions in $\mathcal{A}$ that are invertible form a group $\mathcal{A}^{*}$, which contains an essential subgroup called the set of the multiplicative functions: $\mathcal{M}:=\{f\in\mathcal{A}:f\not\equiv0, f(mn)=f(m)f(n)\text{ for any two relatively prime positive integers }m\text{ and }n\}.$ An interesting result about Dirichlet convolution is the M\"{o}bius inversion formula showing that $f=g*u$ if and only if $f*\mu=g$ where $u\equiv 1$ and $\mu$ is the M\"{o}bius function.

Besides Dirichlet convolution, many other arithmetical convolutions also have long been studied by mathematicians. (See \cite{Subbarao, Lehmer1}.) To prove the results in this paper, we introduce one of them. The $lcm$ convolution of $f$ and $g$ is defined as follows
$$f\circ g=\sum_{[a,b]=n}f(a)g(b).$$
It is obvious that $\mathcal{M}$ is closed under this binary operation. Lehmer \cite{Lehmer1} proved that
\begin{lemma}\label{lem:lehmerconvolution2}
Let $$M(1)=1,\qquad M(n):=\prod_{p^{\alpha}\parallel n}\big((\alpha+1)^{-1}-{\alpha}^{-1}\big).$$
Then $M\circ u=I$ where $I(n)=\lfloor1/n\rfloor$. Thus $f=g\circ u$ if and only if $f\circ M=g.$
\end{lemma}

In another paper, Lehmer \cite{Lehmer2} gave the following formula
\begin{lemma}\label{lem:lehmerconvolution4}
$$(f*u)(g*u)=(f\circ g)*u$$
for any arithmetical functions $f$ and $g$.
\end{lemma}

By the definition of $lcm$ convolution, Lemma \ref{lem:lehmerconvolution2} and Lemma \ref{lem:lehmerconvolution4}, one can easily obtain the following lemmas.

\begin{lemma}
(1) $f\circ\mu=f(1)\mu$ for any $f$ in $\mathcal{A}$.

(2) $f \circ g=I$ if and only if $\Big(\sum_{d\mid n}f(d)\Big)\Big(\sum_{d\mid n}g(d)\Big)=1$ for any positive integer $n$.

(3) $$(f\circ g)(p^{\alpha})=f(p^{\alpha})\sum_{i=0}^{\alpha-1}g(p^i)+g(p^{\alpha})\sum_{j=0}^{\alpha-1}f(p^j)+f(p^{\alpha})g(p^{\alpha})$$
where $p$ is a prime number.

(4) Let $f$ be a multiplicative function. Then $$\sum_{[a,b]=n}\mu(a)f(a)=\prod_{p^2\mid n}(1-f(p))\prod_{p\parallel n}(1-2f(p)).$$ In particular, if $n$ is a perfect square, then $\mu f*u=\mu f\circ u.$
\end{lemma}

\section{proofs of the theorems}

\begin{proof}[Proof of Theorem \ref{th:2.3.1}]
If $\delta=1$, then $a=1$. Thus \eqref{eq:2.3.1} is obvious.

If $\delta=2$, then $m>2$. Assume that $m=2^{\alpha}\prod_{i=1}^{k}p_{i}^{\beta_i}$ where $\{p_i\}$ is a set of prime numbers that is pairwise coprime. By the Chinese remainder theorem,
\begin{eqnarray*}
&\#&\{x\in U_m:x^2\equiv1\pmod m\}\\
=&\#&\{x\in U_{2^{\alpha}}:x^2\equiv1\pmod{2^{\alpha}}\}\prod_{i=1}^{k}\#\{x\in U_{p^{\beta_i}}:x^2\equiv1\pmod{p_{i}^{\beta_i}}\}.
\end{eqnarray*}
Hence,
\begin{equation*}
\#:=\#\{x\in U_m:x^2\equiv1\pmod m\}=\begin{cases}2^k,&\text{ if}\alpha=0\text{ or }1\\2^{k+1},&\text{ if}\ alpha=2\\2^{k+2},&\text{ if}\alpha\geq3.\end{cases}
\end{equation*}

Notice that if $a\not\equiv\pm1\pmod m$, then $ind_m(a)=2$ if and only if $ind_m(-a)=2$. Since $a$ is coprime to $m$, $a\not\equiv -a\pmod m$ if $m>2$. Besides, $ind_m(a)=2$ implies $-a^2\equiv-1\pmod m$. Therefore,
\begin{eqnarray*}
\prod_{\substack{a\in U_m\\ind_m(a)=2}}a&\equiv&(-1)^{\#/2}\pmod m\\
&\equiv&\begin{cases}-1\pmod m,&\text{ if}\alpha=0,1,k=1\text{ or }\alpha=2,k=0,\\1\pmod m,&\text{ otherwise.}\end{cases}
\end{eqnarray*}
The product is congruent to $-1$ modulo $m$ implies $m$ has a primitive root. Thus, \eqref{eq:2.3.1} holds for $\delta=2$.

For $\delta>2$, $ind_m(a)=\delta$ implies $ind_m(a^{-1})=\delta$ where $a^{-1}$ is the multiplicative inverse of $a$ in $U_m$. By the definition of $ind$, $a\not\equiv a^{-1}\pmod m$. Therefore, all integers with $ind$ $delta$ can be arranged in pairs with product $1$ modulo $m$. This proves the theorem.
\end{proof}

\begin{proof}[Proof of Theorem \ref{th:2.3.2}]
(1)is obvious.

By the proof of Theorem \ref{th:2.3.1}, if $a\not\equiv\pm1\pmod m$, then $ind_m(a)=2$ if and only if $ind_m(-a)=2$. Therefore, with the exception of $-1$, all integers with $ind$ $2$ can be arranged in pairs with sum $0$ modulo $m$, which proves (2).

If $4$ divides $\delta$, then $(-a)^{\delta}\equiv1\pmod m$ for $a$ such that $ind_m(a)=\delta$. Let $\delta'$ denote $ind_m(-a)$. If $\delta'<\delta$, then $(-1)^{\delta'}\equiv-1\pmod m$. Otherwise, $a^{\delta'}\equiv1\pmod m$, which contradicts the definition of $ind$. Thus, $\delta'$ is odd and $\delta'\mid\delta$. This implies $a^{2\delta'}\equiv-1\pmod m$. Moreover, $\delta\mid2\delta'$. Hence, $\delta=2\delta'$ which contradicts the assumption that $4$ divides $\delta$. Therefore, $ind_m(a)=\delta$ if and only if $ind_m(-a)=\delta$. For $m>4$, all integers with $ind$ $\delta$ can be arranged in pairs with sum $0$ modulo $m$, which proves (3).
\end{proof}

\begin{proof}[Proof of Theorem \ref{th:2.3.3}]
One can check the congruences in (1), (2) and (3) easily. (4) is a straightforward corollary to Theorem \ref{th:2.3.2}(3).

We now prove the result in (5).

Since $m=p^{\alpha}$, $m$ has a primitive root denoted by $g$. By Lemma \ref{apostolind1}, for any positive integer $d$ that divides $\phi(m)$, there exists integer $a$ such that $ind_m(a)=d$. For any integer $a$ such that $ind_m(a)=\delta$, if $(k,\delta)=1$, then $ind_m(a^k)=\delta$. Notice that $a^i\equiv a^j\pmod m$ if and only if $i\equiv j\pmod{\delta}$. Thus, there exists at least $\phi(\delta)$ integers with $ind$ $\delta$ modulo $m$. Since $n=\sum_{d\mid n}\phi(d)$, there exists exactly $\phi(\delta)$ integers with $ind$ $\delta$ modulo $m$.

Hence, for any integer $a$ such that $ind_m(a)=\delta$,

\begingroup
\allowdisplaybreaks
\begin{eqnarray*}
\sum_{\substack{a\in U_{p^{\alpha}}\\ind_{p^{\alpha}}(a)=\delta}}a&\equiv&\sum_{\substack{j=1\\(j,\delta)=1}}^{\delta}a^j\equiv \sum_{j=1}^{\delta}a^j\sum_{\substack{d\mid j\\d\mid \delta}}\mu(d)\\
&\equiv&\sum_{d\mid\delta}\mu(d)\sum_{\substack{j=1\\d\mid j}}^{\delta}a^j\equiv\sum_{d\mid\delta}\mu(d)\sum_{j=1}^{\delta/d}a^{jd}\\
&\equiv&\sum_{d\mid\delta}\mu(d)a^d\frac{a^{\delta}-1}{a^d-1}\equiv (a^{\delta}-1)\sum_{d\mid\delta}\mu(d)(1+\frac{1}{a^d-1})\\
&\equiv&(a^{\delta}-1)\sum_{d\mid\delta}\mu(d)+\sum_{d\mid\delta}\mu(d)\frac{a^{\delta}-1}{a^d-1}\\
&\equiv&\sum_{d\mid\delta}\mu(d)\frac{a^{\delta}-1}{a^d-1}\equiv\sum_{\substack{d\mid\delta\\a^d\equiv1\pmod p}}\mu(d)\frac{a^{\delta}-1}{a^d-1}\\
&\equiv&\sum_{\substack{d\mid\delta\\g^{\phi(m)d/\delta}\equiv1\pmod p}}\mu(d)\frac{g^{\phi(m)}-1}{g^{\phi(m)d/\delta}-1}\pmod m.
\end{eqnarray*}
\endgroup
Notice that the following four are equivalent
$$g^{\phi(m)d/\delta}\equiv1\pmod p,$$ $$\phi(m)d/\delta\equiv0\pmod{p-1},$$ $$(p-1)d/\delta\equiv0\pmod{p-1},$$ $$(\delta/d,p-1)=1.$$
Thus,
\begin{eqnarray*}
\sum_{\substack{a\in U_{p^{\alpha}}\\ind_{p^{\alpha}}(a)=\delta}}a&\equiv&\sum_{\substack{d\mid\delta\\(\delta/d,p-1)=1}}\mu(d)\sum_{j=1}^{\delta/d}\binom{\delta/d}{j}(g^{\phi(m)d/\delta}-1)^{j-1}\\
&\equiv&\sum_{d\mid p^{ord_p(\delta)}}\mu(d(\delta,p-1))\sum_{j=1}^{\delta/(d(\delta,p-1))}\binom{\delta/(d(\delta,p-1))}{j}(g^{\phi(m)d(\delta,p-1)/\delta}-1)^{j-1}\\
&\equiv&\mu((\delta,p-1))\sum_{d\mid p^{ord_p(\delta)}}\mu(d)\\
&&\times\sum_{j=1}^{\delta/(d(\delta,p-1))}\binom{\delta/(d(\delta,p-1))}{j}(g^{\phi(m)d(\delta,p-1)/\delta}-1)^{j-1}\pmod m.
\end{eqnarray*}
(5) holds apparently for $ord_p(\delta)=0$. For $r:=ord_p(\delta)>0$,
\begin{eqnarray*}
\sum_{\substack{a\in U_{p^{\alpha}}\\ind_{p^{\alpha}}(a)=\delta}}a&\equiv&\mu((\delta,p-1))\{\sum_{j=1}^{p^r}\binom{p^r}{j}(g^{\phi(p^{\alpha-r})}-1)^{j-1}-\sum_{j=1}^{p^{r-1}}\binom{p^{r-1}}{j}(g^{\phi(p^{\alpha-r+1})}-1)^{j-1}\}\\
&\equiv&\mu((\delta,p-1))(p^r-p^{r-1})\equiv\mu((\delta,p-1))\phi(p^r)\pmod m
\end{eqnarray*}
which completes the proof.
\end{proof}

\begin{proof}[Proof of Theorem \ref{th:2.3.4}]
Since $p$ and $q$ are two relatively prime prime numbers, there exists two integers $x$ and $y$ such that $px+qy=1$. Thus,
\begin{eqnarray*}
\sum_{\substack{a\in U_{pq}\\ind_{pq}(a)=\delta}}a&\equiv&\sum_{\substack{a=1\\ind_{p}(a)=\delta_1,ind_{q}(a)=\delta_2\\ [\delta_1,\delta_2]=\delta}}^{pq}a \equiv \sum_{\substack{1\leq i\leq p,1\leq j\leq q\\ind_{p}(i)=\delta_1,ind_{q}(j)=\delta_2\\ [\delta_1,\delta_2]=\delta}}(pxj+qyi)\\
&\equiv&\sum_{[\delta_1,\delta_2]=\delta}\sum_{\substack{i=1\\ind_{p}(i)=\delta_1}}^{p}\sum_{\substack{j=1\\ind_{q}(j)=\delta_2}}^{q}i\equiv\sum_{[\delta_1,\delta_2]=\delta}\sum_{\substack{i=1\\ind_{p}(i)=\delta_1}}^{p}i\sum_{\substack{j=1\\ind_{q}(i)=\delta_2}}^{p}1\\
&\equiv&\sum_{\substack{[\delta_1,\delta_2]=\delta\\ \delta_1\mid p-1,\delta_2\mid q-1}}\mu(\delta_1)\phi(\delta_2)\pmod p.
\end{eqnarray*}

Let $$\chi_{n}(a):=\begin{cases}1,&\text{ if }a\mid n,\\0,&\text{ if }a\nmid n,\\\end{cases}.$$
Then
\begin{eqnarray*}
\sum_{\substack{a\in U_{pq}\\ind_{pq}(a)=\delta}}a&\equiv&\sum_{[\delta_1,\delta_2]=\delta}\mu(\delta_1)\chi_{p-1}(\delta_1)\phi(\delta_2)\chi_{q-1}(\delta_2)\\
&\equiv&\mu\chi_{p-1}\circ\phi\chi_{q-1}(\delta)\\
&\equiv&\sum_{d\mid \delta}\mu(\delta/d)\sum_{x\mid d}\mu(x)\chi_{p-1}(x)\sum_{y\mid d}\phi(y)\chi_{p-1}(y)\\
&\equiv&\sum_{\delta_1\mid(\delta,p-1)}\mu(\delta_1)\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=\delta}}\phi(\delta_2)\pmod p.
\end{eqnarray*}

If $(\frac{\delta}{(\delta,p-1)},p-1)\neq1$, then there exists some prime $r$ that divides $(\frac{\delta}{(\delta,p-1)},p-1)$.
Hence, $ord_r(\delta)>ord_r((\delta,p-1))$. Moreover, $ord_r(\delta)>ord_{r}(p-1)>0.$

Let $\delta_1$ be a divisor of $(\delta,p-1)/r$. Then $ord_{r}(\delta_1)<ord_{r}(p-1)<ord_{r}(\delta)$. Thus, $ord_{r}(r\delta_1)<ord_{r}(\delta).$ Therefore, $[\delta_1,\delta_2]=\delta$ if and only if $[r\delta_1,\delta_2]=\delta$, which leads to the following
$$\sum_{\delta_1\mid(\delta,p-1)}\mu(\delta_1)\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=\delta}}\phi(\delta_2)=\sum_{\delta_1\mid(\delta,p-1)/r}(\mu(\delta_1)+\mu(r\delta_1))\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=\delta}}\phi(\delta_2)=0.$$
Therefore we proved the case of $(\frac{\delta}{(\delta,p-1)},p-1)>1$.

Assume $(\frac{\delta}{(\delta,p-1)},p-1)=1$. If $(\delta,p-1)=1$, then
$$\sum_{\delta_1\mid(\delta,p-1)}\mu(\delta_1)\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=\delta}}\phi(\delta_2)=\phi(\delta).$$
This proves the theorem.

If $(\delta,p-1)>1$, then $ord_{r}(\delta)\leq ord_{r}(p-1)$ for any prime $r$ that divides $(\delta,p-1)$. Thus, $(\frac{\delta}{(\delta,p-1)},\delta_1)=1$. Hence,
\begin{eqnarray*}
\sum_{\delta_1\mid(\delta,p-1)}\mu(\delta_1)\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=\delta}}\phi(\delta_2)&=&\sum_{\delta_1\mid(\delta,p-1)}\mu(\delta_1)\sum_{\substack{\delta_2\mid q-1\\ [\delta_1,\delta_2]=(\delta,p-1)}}\phi(\delta_2)\phi\Big(\frac{\delta}{(\delta,p-1)}\Big)\\
&=&\frac{\phi(\delta)}{\phi((\delta,p-1))}\mu\circ\phi\chi_{q-1}((\delta,p-1))\\
&=&\frac{\phi(\delta)}{\phi((\delta,p-1))}\mu((\delta,p-1))\phi\chi_{q-1}(1)\\
&=&\frac{\phi(\delta)}{\phi((\delta,p-1))}\mu((\delta,p-1)).
\end{eqnarray*}
This completes the proof.
\end{proof}

\begin{proof}[Proof of Theorem \ref{th:2.3.5}]

By Theorem \ref{th:2.3.3} (5) and the methods mentioned in the proof of Theorem \ref{th:2.3.4}, we can obtain the following congruence
\begin{equation}\label{final}
\sum_{\substack{a\in U_{m}\\ind_{m}(a)=\delta}}a\equiv \big(\mu\chi_{p-1}\circ\phi\chi_{\phi(p_1^{\alpha_1})}\circ\cdots\circ\phi\chi_{\phi(p_k^{\alpha_k})}\big)(\delta)\pmod{p^{\alpha}}.
\end{equation}
Now, we will show how \eqref{final} leads to \eqref{eq:th:2.3.5.2}. Let $rad(n):=\prod_{p\mid n}p$. Then for any positive integer $a$ and multiplicative function $g$,

\begin{eqnarray*}
\mu\chi_{a}\circ g(n)&=&\prod_{p^{\alpha}\parallel n}\{I((p^{\alpha},a))g*u(p^{\alpha})-I((p^{\alpha-1},a))g*u(p^{\alpha-1})\}\\
&=&I((\frac{n}{rad(n)},a))\prod_{\substack{p\parallel n\\p\mid a}}(-1)g*u(1)\prod_{\substack{p^{\alpha}\parallel n\\p\nmid a}}(g*u(p^{\alpha})-g*u(p^{\alpha-1}))\\
&=&I((\frac{n}{rad(n)},a))\prod_{\substack{p\parallel n\\p\mid a}}(-1)\prod_{p^{\alpha}\parallel \frac{n}{(n,a)}}(\mu(1)g*u(p^{\alpha})+\mu(p)g*u(p^{\alpha-1}))\\
&=&\mu((n,a))I((\frac{n}{(n,a)},a))\prod_{p^{\alpha}\parallel \frac{n}{(n,a)}}\mu*g*u(p^{\alpha})\\
&=&\mu((n,a))I((\frac{n}{(n,a)},a))g(\frac{n}{(n,a)})
\end{eqnarray*}

For any positive integers $b_1,$ $\cdots$, $b_k$,
\begin{eqnarray*}
\phi\chi_{b_1}\circ\cdots\circ\phi\chi_{b_k}(n)&=&\prod_{p^{\alpha}\parallel n}\{\sum_{j=0}^{\alpha}\mu(p^{\alpha-j})\prod_{i=1}^{k}(\phi\chi_{b_i}*u)(p^j)\}\\
&=&\prod_{p^{\alpha}\parallel n}\{\sum_{j=0}^{\alpha}\mu(p^{\alpha-j})\prod_{i=1}^{k}(\phi\chi_{b_i}*u)(p^j)\}\\
&=&\prod_{p^{\alpha}\parallel n}\{\prod_{i=1}^{k}(\phi\chi_{b_i}*u)(p^{\alpha})-\prod_{i=1}^{k}(\phi\chi_{b_i}*u)(p^{\alpha-1})\}\\
&=&\prod_{p^{\alpha}\parallel n}\prod_{ord_p(b_i)<\alpha}p^{ord_p(b_i)}\prod_{ord_p(b_i)\geq\alpha}(p^{\alpha}-p^{\alpha-1}).
\end{eqnarray*}

Combining these two formulas, we have
\begin{eqnarray*}
&&\big(\mu\chi_{p-1}\circ\phi\chi_{\phi(p_1^{\alpha_1})}\circ\cdots\circ\phi\chi_{\phi(p_k^{\alpha_k})}\big)(\delta)\\
&=&\mu((\delta,p-1))I((\frac{\delta}{(\delta,p-1)},p-1))F(\phi(p_1^{\alpha_1}),\cdots,\phi(p_k^{\alpha_k});\frac{\delta}{(\delta,p-1)}).
\end{eqnarray*}
\end{proof}

\subsection*{Acknowledgments}
This work is supported by the National Natural Science Foundation of China (Grant No. 11571303) and the Natural Science Foundation of Zhejiang Province (No. LY18A010016).

\bibliographystyle{amsplain}