\section{Guided Learning through Successive Functional Gradient Optimization} 
\label{sec:theory}

In this section, after presenting the framework in general terms, 
we develop concrete algorithms and analyze them.  

\subsection{Framework} 

We first describe the framework in general terms so that the models to be 
trained are not limited to parameterized ones. 
Let $f$ be the model we are training.  
Starting from some initial $f$, training proceeds by repeating the following: 
\begin{enumerate}[itemsep=1pt,topsep=1pt]
  \item Generate a {\em guide function} $f^*$ by applying 
        functional gradient descent for reducing the loss to the current model $f$, 
        so that $f^*$ is an improvement over $f$ in terms of loss 
        but not too far from $f$.  
  \item Move the model $f$ in the direction of the guide function $f^*$
        according to some distance measure. 
\end{enumerate}

We use the Bregman divergence $\Bdiv$, defined in \eqref{eq:bdiv}, 
for representing the distances between models.

\tightpara{Step 1: Guide going ahead}
We formulate Step 1 %
as 
\begingroup
\thinmuskip=0mu
\medmuskip=0mu
\thickmuskip=0mu
\begin{align}
    f^*(x,y) &:= \arg\min_q \left[ 
        \Bdiv(q,f(x)) + \alpha \nabla \loss_y(f(x))^\top q
    \right ] , \label{eq:mirror0}  
\end{align}
\endgroup
where $\alpha$ is a meta-parameter.  
The second term pushes the guide function towards the direction of reducing loss, 
and the first term pulls back the guide function towards the current model $f$.  
Thus, $f^*$ is ahead of $f$ but not too far ahead.  
Note that we use the knowledge of the true output $y$ here; therefore, $f^*$ takes 
$y$ as the second argument.  
The function value for each data point $(x,y)$ can be found approximately 
by solving the optimization problem by SGD if there is no analytical solution. 
Also, this formulation is equivalent to finding $f^*$ such that 
\begin{align}
  \nabla \Bref(f^*(x,y)) = \nabla \Bref(f(x)) - \alpha \nabla \loss_y(f(x)) . 
  \label{eq:mirror1} 
\end{align}
This is mirror descent (see e.g., \cite{convex15}) performed in a function space.  

Due to the relation of the Bregman divergence to the Hessian matrix stated in \eqref{eq:bdivApprox}, 
\eqref{eq:mirror0} implies that 
\begin{align}
  f^*(x,y) \approx f(x) - \alpha (\hessian{\Bref(f(x))})^{-1} \nabla \loss_y(f(x)) .
  \label{eq:step}
\end{align}

Therefore, if we set $\Bref(f)=\loss_y(f)$, \eqref{eq:step} becomes 
\begin{align}
  f^*(x,y) \approx f(x) - \alpha (\hessian{\loss_y(f(x))})^{-1} \nabla \loss_y(f(x)), 
  \label{eq:newton}
\end{align}
which is approximately a second-order functional gradient step 
(one step of the relaxed Newton method) with step-size $\alpha$ for minimizing the loss.  

If we set $\Bref(f)=\frac{1}{2}\|f\|^2$, then the optimization problem \eqref{eq:mirror0} 
has an analytical solution
\[
  f^*(x,y) = f(x) - \alpha \nabla \loss_y(f(x)), 
\]
which is a first-order functional gradient step with step-size $\alpha$
for minimizing the loss.  

\tightpara{Taking $m$ steps in Step 1}
For further generality, let us also consider $m$ steps of functional gradient descent by 
extending $f^*$ in \eqref{eq:mirror0} to $f^*_m$ recursively defined as follows. 
\begin{align*}
    f^*_0(x,y) &:= f(x) \\
    f^*_{i+1}(x,y) &:= \arg\min_q \left[ 
        \Bdiv(q,f^*_i(x)) + \alpha \nabla \loss_y(f^*_i(x))^\top q
    \right ]. %
\end{align*}
Then, 
in parallel to \eqref{eq:step}, we have 
\begin{align*}
  f^*_{i+1}(x,y) \approx f^*_i(x) - \alpha (\hessian{\Bref(f^*_i(x))})^{-1} \nabla \loss_y(f^*_i(x)) .
\end{align*}

\tightpara{Step 2: Following the guide}  
Using the Bregman divergence $\Bdiv$, we formulate Step 2 above as an update of the model $f$ 
to reduce 
\begin{align}
    \meanOverS{ \Bdiv(f(x), f^*(x,y)) } + R(f)
    \label{eq:fit}
\end{align}
so that the model $f$ approaches the guide function $f^*$ 
in terms of the Bregman divergence.  
$R(f)$ is a regularization term.  

\tightpara{Parameterization}
Although there can be many variations of this scheme, in this work, 
we parameterize the model $f$ so that we can train neural networks.  
Thus, we replace $f(x)$ by $\net{\theta}{x}$ with parameter $\theta$.  
This does not affect Step 1, and to reduce \eqref{eq:fit} in Step 2, 
we repeatedly update the model parameter $\theta$ by descending the stochastic gradient 
\begin{align}
\nabla_\theta \left[ 
   \meanOverB { \Bdiv(\net{\theta}{x}, f^*(x,y)) } + R(\theta)
\right] , 
\end{align}
where $B$ is a mini-batch sampled from a training set $S$.  

\subsection{Algorithms}

Putting everything together, we obtain Algorithm \ref{alg:Algm}, 
which performs mirror descent in a function space 
in Line 3.  
We call it (and its derivatives) a method of 
{\em GUided Learning through successive Functional gradient optimization (GULF)}.   
We now instantiate function $\Bref$ used by the Bregman divergence $\Bdiv$ 
to derive concrete algorithms.  In general we allow $\Bref$ to vary
for each data point. That is, it may depend on $(x,y)$. 
Here we use two functions discussed above, which correspond to 
the first-order and the second-order methods, respectively; 
however, note that choice of $\Bref$ is not limited to these two.  

\begin{algorithm}[h]
\begin{small}
   \caption{
      GULF\ in the most general form. 
      {\bfseries Input:} %
       $\theta_0$, 
       training set $S$. %
       Meta-parameters: $m$, $\alpha$, $ T $. 
      {\bfseries Output:} $\theta_ T $. 
   }
   \label{alg:Algm}
\begin{algorithmic}[1]
   \STATE $\theta  \leftarrow  \theta_0$
   \FOR{$t=0$ {\bfseries to} $ T -1$}  %
      \STATE Define             
               $\tarffx{m}$ by: $\tarffx{0}(x,y) := \net{\theta_t}{x}$,
               $\tarffx{i+1}(x,y) := $\\ $~~~~~~~~~~~~\arg\min_q \left[ 
                  \Bdiv(q,\tarffx{i}(x,y)) + \alpha \nabla \loss_y(\tarffx{i}(x,y))^\top q
                \right ]$        
      \REPEAT 
         \STATE Sample a mini-batch $B$ from $S$.
         \STATE  Update $\theta$ by descending the stochastic gradient \\
             $\nabla_\theta \left[ 
               \meanOverB{ \Bdiv(\net{\theta}{x}, \tarffx{m}(x,y)) } 
               + R(\theta)
             \right]$ \\
             $~~~~$for optimizing \\
             $Q_t(\theta) :=     \meanOverS{ \Bdiv(\net{\theta}{x}, \tarffx{m}(x,y)) } 
               + R(\theta)$.
      \UNTIL{some criteria are met}
      \STATE $\theta_{t+1}  \leftarrow  \theta$
   \ENDFOR
\end{algorithmic}
\end{small}
\end{algorithm}

\tightpara{\ourssq\ (1st-order, Algorithm \ref{alg:Algsq})} 
With $\Bref(u)=\frac{1}{2}\| u \|^2$, %
we obtain Algorithm \ref{alg:Algsq}.  Derivation is straightforward.  
This algorithm performs $m$ steps of first-order functional gradient descent (Line 3)
to push the guide function ahead of the current model 
and then let the model follow the guide by reducing the 2-norm between them.   

\begin{algorithm}[h]
\begin{small} 
   \caption{\ourssq\ ($\Bref(u)=\frac{1}{2}\|u\|^2$): 
       {\bfseries Input:} %
       $\theta_0$, 
       training set $S$. %
       Meta-parameters: $m$, $\alpha$, $ T $.    
       {\bfseries Output:} $\theta_ T $. 
   }
   \label{alg:Algsq}
\begin{algorithmic}[1]
   \STATE $\theta  \leftarrow  \theta_0$
   \FOR{$t=0$ {\bfseries to} $ T -1$}   
      \STATE Define $\tarffx{m}$ by: $\tarffx{0}(x,y) = \net{\theta_t}{x}$, \\
      $~~~~~~~~~~~~~~~~~~~~~~~~$$\tarffx{i+1}(x,y) = \tarffx{i}(x,y) - \alpha \nabla \loss_y(\tarffx{i}(x,y))$   
      \REPEAT
         \STATE Sample a mini-batch $B$ from $S$. 
         \STATE Update $\theta$ by descending the stochastic gradient \\
             $~~~~~~\nabla_\theta \left[ 
               \meanOverB{ \frac{1}{2}\| \net{\theta}{x} - \tarffx{m}(x,y) \|^2 } 
               + R(\theta)
             \right]$ 
      \UNTIL{some criteria are met}
      \STATE $\theta_{t+1}  \leftarrow  \theta $
   \ENDFOR
\end{algorithmic}
\end{small}
\end{algorithm}


\begin{algorithm}[h]
\begin{small}
   \caption{\oursL\ ($\Bref(p)=\loss_y(p)$): 
       {\bfseries Input:} %
       $\theta_0$, 
       training set $S$. %
       Meta-parameters: $\alpha \in (0,1)$, $ T $.   
       {\bfseries Output:} $\theta_ T $. 
       Notation: $\netprm=\net{\prm}{x}$ and $\netprmt=\net{\prm_t}{x}$.
   }
   \label{alg:AlgL}
\begin{algorithmic}
   \STATE $\theta  \leftarrow  \theta_0$
   \FOR{$t=0$ {\bfseries to} $ T -1$}  
      \REPEAT 
         \STATE Sample a mini-batch $B$ from $S$. 
         \STATE Update $\theta$ by descending the stochastic gradient \\
             $\nabla_\theta \left[ 
               \meanOverB{ \Bdivx{\loss_y}(\netprm, \netprmt) + \alpha \nabla \loss_y(\netprmt)^\top \netprm }               
               + R(\theta)
             \right]$ 
      \UNTIL{some criteria are met}
      \STATE $\theta_{t+1}  \leftarrow  \theta$
   \ENDFOR
\end{algorithmic}
\end{small}
\end{algorithm}

\tightpara{\oursL\ (2nd order, Algorithm \ref{alg:AlgL})} 

We consider the case of $\Bref(p)=\loss_y(p)$ (i.e., $\Bref$ returns loss given prediction $p$).  
\eqref{eq:newton} has shown that in this case 
Step 1 becomes approximately the second-order functional gradient descent.  
Also, with this choice of $\Bref$, 
Algorithm \ref{alg:Algm}\ can be converted to a simpler form where 
we do not have to compute the values of the guide function $\tarffx{m}$ explicitly, 
and where we have one fewer meta-parameter.  
This simpler form is shown in Algorithm \ref{alg:AlgL}, which has the following 
relationship to Algorithm \ref{alg:Algm}.  

\begin{proposition}
\label{pr:m}
When $\Bref(p) = \loss_y(p)$ that returns loss given prediction $p$, 
Algorithm \ref{alg:Algm}\ with $\alpha=\gamma$ is equivalent to
Algorithm \ref{alg:AlgL}\ with $\alpha=1 - (1-\gamma)^m$. 
\end{proposition}
The proofs are all provided in the supplementary material. 

To simplify notation, let $\netprm=\net{\theta}{x}$, which is the model that we are updating, 
and $\netprmt=\net{\theta_t}{x}$, which is a model that was frozen when time changed from $t-1$ to $t$.  
In the stage associated with time $t$, Algorithm \ref{alg:AlgL}\ minimizes
\begin{equation}
  \meanOverS{ \Bdivx{\loss_y}(\netprm,\netprmt) + \alpha \nabla \loss_y(\netprmt)^\top \netprm} + R(\theta)
  \label{eq:algLobj0}
\end{equation}
approximately through mini-batch SGD. 
The second term $\alpha \nabla \loss_y(\netprmt)^\top \netprm$ pushes the model $\netprm$ 
towards the direction of reducing loss, 
and the first term $\Bdivx{\loss_y}(\netprm,\netprmt)$ pulls it back towards the frozen model $\netprmt$.  
With a certain family  of loss functions, \eqref{eq:algLobj0} can be further transformed 
as follows. 
\begin{proposition}
\label{pr:distill}
Let $y$ be a vector representation 
such as a $K$-dim vector representing $K$ classes.  
Assume that the gradient of the loss function can be expressed as 
\begin{equation}
  \nabla L(f,y) = \nabla \loss_y(f) = p(f) - y  \label{eq:lossCond}
\end{equation}
with $p(f)$ not depending on $y$. Let
\begin{align}
J_t(\theta) =& \meanOverS{\Bdivx{\loss_y}(\netprm, \netprmt) + \alpha \nabla \loss_y(\netprmt)^\top
\netprm} 
\label{eq:algLobj} \\
J'_t(\theta) =& \meanOverS{(1-\alpha)L(\netprm, p(\netprmt)) + \alpha
\loss_y(\netprm)} 
\label{eq:distillobj}              
\end{align}
Then we have
\[
J_t(\theta) = J'_t(\theta) + c_t~,
\]
where $c_t$ is independent of $\theta$.
This implies that
\[
\arg\min_\theta \left[ J_t(\theta) + R(\theta) \right]
=
\arg\min_\theta \left[ J'_t(\theta) + R(\theta) \right] .
\]
\end{proposition}
Both the cross-entropy loss and squared loss satisfy \eqref{eq:lossCond}.  
In particular, when $\loss_y(f)$ is the cross-entropy loss, 
$p(f)$ becomes the soft max function.  
In this case, \eqref{eq:distillobj} is the {\em distillation} formula 
with the frozen model $\netprmt$ playing the role of a cumbersome source model, 
and therefore, the parameter update rule of Algorithm \ref{alg:AlgL}\ involving \eqref{eq:algLobj} 
becomes that of distillation.  
Thus, Algorithm~\ref{alg:AlgL} can be regarded as a generalization of
self-distillation for arbitrary loss functions.

\subsection{Convergence Analysis}

Let us define {\em $\alpha$-regularized loss} 
\begin{equation}
  \reglossx{\alpha}(\theta) := \meanOverSp{ \myloss{\net{\theta}{x}} } + \frac{1}{\alpha}R(\theta).  \label{eq:reglossa}
\end{equation}
The following theorem shows that 
Algorithm \ref{alg:Algm} with step-size $\alpha$ always approximately reduces the $\alpha$-regularized loss 
if $\alpha$ is appropriately set.  

\newcommand{\tilde{\theta}}{\tilde{\theta}}
\begin{theorem}
\label{thm:L-conv}
In the setting of Algorithm \ref{alg:Algm}\ with $m=1$, assume that there exists $\beta>0$
such that $\Bdivx{h}(f,f') \geq \beta \Bdivx{\loss_y}(f,f')$ for any $f$ and $f'$, and 
assume that $\alpha \in (0,\beta]$.  
Assume also that $Q_t(\theta)$ defined in Algorithm \ref{alg:Algm}\ is $1/\eta$ smooth in $\theta$:
$
\|\nabla Q_t(\theta) - \nabla Q_t(\theta')\| \leq (1/\eta) \|\theta -
\theta'\| .
$

Assume that $\theta_{t+1}$ is an improvement of $\theta_t$ with respect to minimizing $Q_t$ so that 
$Q_t(\theta_{t+1}) \leq Q_t(\tilde{\theta})$ ,  
where
\begin{align}
\tilde{\theta} = \theta_t - \eta \nabla Q_t(\theta_t) . 
\label{eq:algLupd}
\end{align}
Then we have 
\[
  \reglossx{\alpha}(\theta_{t+1}) \le \reglossx{\alpha}(\theta_t) - \frac{\alpha \eta}{2}
  \|\nabla \reglossx{\alpha}(\theta_{t})\|^2 .
\]
\end{theorem}
For Algorithm \ref{alg:AlgL}, we have $h(\cdot)=\loss_y(\cdot)$ and thus $\beta=1$, 
leading to $\alpha\in (0,1]$.  
\eqref{eq:algLupd} is the parameter update step of Algorithm \ref{alg:Algm}\ 
except that the algorithm stochastically estimates the mean over $S$ 
from a mini-batch $B$ sampled from $S$.  
Therefore, 
the theorem indicates that each stage (corresponding to $t$) of the algorithm 
approximately reduces the $\alpha$-regularized loss $\reglossx{\alpha}(\theta)$.  
In other words, while the guide function changes from stage to stage, 
a quantity that does {\em not} depend on the guide function 
goes down throughout training, namely, the $\alpha$-regularized loss $\reglossx{\alpha}$.  

Furthermore, we obtain from Theorem~\ref{thm:L-conv} that
\[
\frac{1}{T}
\sum_{t=0}^{T-1}   \|\nabla \reglossx{\alpha}(\theta_{t})\|^2 
\leq \frac{2( \reglossx{\alpha}(\theta_0)-\reglossx{\alpha}(\theta_{T}))}{\alpha \eta T} .
\]
Assuming $\reglossx{\alpha}(\theta) \ge 0$, this implies that as $ T $ goes to infinity, 
the right-hand side goes to zero, and so Algorithm \ref{alg:AlgL}\ converges with 
$\nabla \reglossx{\alpha}(\theta_{ T }) \to 0$.
Therefore, when $ T $ is sufficiently large, $\theta_{ T }$ finds a
stationary point of $\reglossx{\alpha}$. 

The convergence result indicates that having a regularization term $R(\theta)$ in the algorithm 
effectively causes minimization of the $\alpha$-regularized loss.  
However, our empirical results (shown later) indicate that 
GULF\ models are very different from standard models trained directly to
minimize the $\alpha$-regularized loss.  
For example, standard models trained with $\reglossx{0.01}$ suffers from severe 
underfitting, but GULF\ model with $\alpha$=0.01 produces high performance.  
This is because each
step of guided learning tries to find a good solution which is near
the previous solution (guidance). The complexity of each iterate is better controlled,
and hence this approach leads to better generalization performance. 
We will come back to this point in the next section.  



\section{On the Empirical Study} 

In this section, we first 
provide experimental details and additional figures regarding the experiments 
reported in the main paper, 
and then we report additional experiments using text data.  
Our code is provided at a repository under {\tt github.com/riejohnson}.  

\subsection{Details of the experiments in the main paper}

\subsubsection{CIFAR10, CIFAR100, and SVHN}

This section describes the experimental details of all but the ImageNet experiments.  

The mini-batch size was set to 128.  We used momentum 0.9.  
The following learning rate scheduling was used: 
200K steps with $\eta$, 40K steps with $0.1\eta$, and 40K steps with $0.01\eta$.  
The initial learning rate $\eta$ was set to 0.1 on CIFAR10/100 
and 0.01 on SVHN, following \cite{Wresnet16}. 
The weight decay $\lambda$ was 0.0001 except that 
it was 0.0005 for (CIFAR100, WRN-28-10) and SVHN.  

We used the standard mean/std normalization on all 
and the standard shift and horizontal flip image augmentation on CIFAR10/100.  

We report the median of three runs with three random seeds.  
The meta-parameters were chosen based on the performance on the development set.  
All the results were obtained by using 
only the `train' portion (shown in Table \ref{tab:data-stat} of the main paper) 
of the official training set as training data.  

For label smoothing, the amount of probability taken away from the true class was chosen 
from $\{0.1, 0.2, 0.3, 0.4\}$.  

To obtain the results reported in Table \ref{tab:c10-c100-svhn} (with smaller networks), 
$ T $ was fixed to 25 for CIFAR10/100, and 15 for SVHN.  
$\alpha$ for ini:random\ was fixed to 0.3. 
For ini:base, we chose $\alpha$ from $\{ 0.3, 0.01 \}$.  
We excluded $\alpha=0.01$ for ini:random, as it takes too long. 
When dropout was applied in the SVHN experiments, 
the dropout rate was set to 0.4, following \cite{Wresnet16}. 
To obtain the results reported in Table \ref{tab:c10-c100-k10} (with larger networks), 
$ T $ was fixed to 1.  
For \oursL, $\alpha$ was chosen from $\{ 0.3, 0.01 \}$. 
For \ourssq, $\alpha$ was fixed to 0.3, and $m$ 
(the number of functional gradient steps) was chosen from 
$\{ 1, 2, 5 \}$.  
On CIFAR datasets, the choice of $\alpha$ or $m$ did not make much difference, 
and the chosen values tended to vary among the random seeds.  
On SVHN, $\alpha$=0.01 tended to be better when no dropout was used, 
and 0.3 was better when dropout was used.  

To perform random initialization of the parameter for ini:random\ and 
the baseline methods, we used  
Kaiming normal initialization \cite{HZRS15}, following the previous work. 

\subsubsection{ImageNet}

Each stage of the training for ImageNet followed the code used for training the 
pre-trained models provided as part of TorchVision: 
{\tt https://github.com/pytorch/examples/blob/master/imagenet/main.py}.  
That is, for both ResNet-50 and WRN-50-2, 
the learning rate was set to 
$\eta$, $0.1\eta$, and $0.01\eta$ for 30 epochs each, 
i.e., 90 epochs in total, and the initial rate $\eta$ was set to 0.1. 
The mini-batch size was set to 256, and the weight decay was set to 0.0001. 
The momentum was 0.9.  
$\alpha$ was fixed to 0.5. 
We used two GPUs for ResNet-50 and four GPUs for WRN-50-2. 

We used the standard mean/std normalization 
and the standard image augmentation for ImageNet -- 
random resizing, cropping and horizontal flip, 
which is the same data augmentation scheme 
as used for training the pre-trained models provided as part of TorchVision. 

\clearpage
\subsection{Additional figures}

Figure \ref{fig:c100k1-joint-terr-with-legend} shows test error (\%) 
in relation to training loss with a small ResNet on CIFAR100. 
Additional examples of test-loss curves are shown in Figure \ref{fig:tloss}. 
Figure \ref{fig:c100k1-l2-alpha} shows the parameter size $\|\theta_t\|^2$ 
in relation to training loss, in the settings of Figure \ref{fig:c100k1-tloss-alpha} 
in the main paper.  

\begin{figure}[ht]
\begin{center}
\centerline{\includegraphics[width=0.25\linewidth]{aug-c100k1-joint-terr-with-legend}}
\vskip -0.15in
\caption{ \label{fig:c100k1-joint-terr-with-legend}
  Test error (\%) in relation to training loss. 
  The arrows indicate the direction of time flow.   
  \oursL. CIFAR100.  ResNet-28.
} 
\end{center}
\vskip 0.2in
\begin{center}
\begin{subfigure}[b]{0.25\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c10k1-tloss-zoom-with-legend}}
\label{fig:c10k1-tloss-zoom}
\vskip -0.1in
\caption{CIFAR10, ResNet-28}
\end{center}
\end{subfigure}%
\quad
\begin{subfigure}[b]{0.25\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{svhn-do-tloss-zoom-with-legend}}
\label{fig:svhn-tloss-zoom}
\vskip -0.1in
\caption{SVHN, WRN-16-4.}
\end{center}
\end{subfigure}%
\quad
\begin{subfigure}[b]{0.25\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c10dense-tloss-with-legend}}
\label{fig:c10dense-tloss}
\vskip -0.1in
\caption{CIFAR10,DenseNetBC-40-12.}
\end{center}
\end{subfigure}%
\end{center}
\vskip -0.35in
\caption{ \label{fig:tloss}
  Additional examples of test loss curves of \oursL.  
  The arrows indicate the direction of time flow.
}
\vskip 0.4in
\centering
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-a1}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-a1}
base-loop\ ($\alpha$=1)
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-a09}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-a09}
  $\alpha$=0.9
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-a03}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-a03}
  $\alpha$=0.3
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-a01}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-a01}
  $\alpha$=0.1
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-a001}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-a001}
  $\alpha$=0.01
}
\end{center}
\end{subfigure}%
\vskip -0.25in
\caption{ \label{fig:c100k1-l2-alpha}
Parameter size $\|\theta_t\|^2$ of ini:base (`$\triangle$') and ini:random(`$\circ$').  
with five values of $\alpha$ (becoming smaller from left to right), 
in relation to training loss. \oursL.  $ T $=25. 
CIFAR100.  ResNet-28. 
Matching figures with Figure \ref{fig:c100k1-tloss-alpha}.  
As $\alpha$ becomes smaller, 
the (potential) meeting point shifts 
further away from the base model. 
The left-most figure is base-loop, which is equivalent to 
$\alpha$=1. 
The arrows indicate the direction of time flow.  
}
\vskip -0.1in
\end{figure}

\input{breg-exp-text}

\section{Proofs}

In the proofs, 
we use abbreviated notation by dropping $x$ and $y$ and making $\theta$ a subscript, e.g., 
we write $\netprm$ for $\net{\prm}{x}$.  

\subsection{Proof of Proposition \ref{pr:m}}

\paragraph{Proposition \ref{pr:m}}
{\em 
  When $\Bref(p) = \loss_y(p)$ that returns loss given prediction $p$, 
  Algorithm \ref{alg:Algm}\ with $\alpha=\gamma$ is equivalent to
  Algorithm \ref{alg:AlgL}\ with $\alpha=1 - (1-\gamma)^m$. 
}
\begin{proof}
From Algorithm \ref{alg:Algm}\ with $\alpha=\gamma$, we have
\begin{align}
   \tarffx{i} = \arg\min_q \left[ 
                \Bdiv(q,\tarffx{i-1}) + \gamma \nabla \loss_y(\tarffx{i-1})^\top q
   \right ] . \label{eq:tarffxi}
\end{align}
From $\Bref(\cdot)=\loss_y(\cdot)$ and \eqref{eq:tarffxi}, we obtain 
\begin{align*}
\nabla\loss_y(\tarffx{i}) &= \nabla\loss_y(\tarffx{i-1}) - \gamma\nabla\loss_y(\tarffx{i-1}) = (1-\gamma)\nabla\loss_y(\tarffx{i-1})
~~\mbox{ for }i=1,\cdots,m .
\end{align*}
Since $\tarffx{0}=\netprmt$, we have 
\begin{align*}
\nabla\loss_y(\tarffx{m}) &= (1-\gamma)^m \nabla\loss_y(\tarffx{0})
= (1-\gamma)^m \nabla\loss_y(\netprmt) , 
\end{align*}
which implies %
\begin{align*}
\nabla_{\netprm} \left[ \Bdiv(\netprm,\tarffx{m}) \right ] 
&=  \nabla\loss_y(\netprm) - \nabla\loss_y(\tarffx{m})
 =  \nabla\loss_y(\netprm) - (1-\gamma)^m \nabla\loss_y(\netprmt) \\
&=  \nabla_{\netprm} \left[ \Bdivx{\loss_y}(\netprm, \netprmt) + ( 1 - (1-\gamma)^m) \nabla \loss_y(\netprmt)^\top \netprm \right]
\end{align*}
and therefore, 
$
\nabla_\theta \left[ \Bdiv(\netprm,\tarffx{m}) \right ] = 
\nabla_\theta \left[ \Bdivx{\loss_y}(\netprm, \netprmt) + ( 1 - (1-\gamma)^m) \nabla \loss_y(\netprmt)^\top \netprm \right]
$.  The rest is trivial. 
\end{proof}

\subsection{Proof of Proposition \ref{pr:distill}}
\paragraph{Proposition \ref{pr:distill}}
{\em 
Let $y$ be a vector representation 
such as a $K$-dim vector representing $K$ classes.  
Assume that the gradient of the loss function can be expressed as 
\[
  \nabla L(f,y) = \nabla \loss_y(f) = p(f) - y  
\]
with $p(f)$ not depending on $y$. Let
\begin{align*}
J_t(\theta) =& \meanOverS{\Bdivx{\loss_y}(\netprm, \netprmt) + \alpha \nabla \loss_y(\netprmt)^\top
\netprm} \\
J'_t(\theta) =& \meanOverS{(1-\alpha)L(\netprm, p(\netprmt)) + \alpha
\loss_y(\netprm)} 
\end{align*}
Then we have
\[
J_t(\theta) = J'_t(\theta) + c_t,
\]
where $c_t$ is independent of $\theta$.
This implies that
\[
\arg\min_\theta \left[ J_t(\theta) + R(\theta) \right]
=
\arg\min_\theta \left[ J'_t(\theta) + R(\theta) \right] .
\]
} %
\begin{proof}
\begin{align*}
\nabla_{\netprm} \left[ \Bdivx{\loss_y}(\netprm, \netprmt) + \alpha \nabla \loss_y(\netprmt)^\top \netprm \right ]
&= \nabla\loss_y(\netprm) - (1-\alpha)\nabla\loss_y(\netprmt) \\
&= (p(\netprm)-y) - (1-\alpha)(p(\netprmt)-y) \\
&= (1-\alpha)(p(\netprm) - p(\netprmt)) + \alpha( p(\netprm) - y ) \\
&= \nabla_{\netprm} \left[ (1-\alpha)L(\netprm, p(\netprmt)) +
  \alpha\loss_y(\netprm) \right ] .
\end{align*}  
This implies that $\nabla J_t(\theta) = \nabla
J'_t(\theta)$. Therefore $J_t(\theta) - J'_t(\theta)$ is independent of $\theta$. 
\end{proof}

\subsection{Proof of Theorem \ref{thm:L-conv}}

\paragraph{Theorem \ref{thm:L-conv}} 
{\em 
In the setting of Algorithm \ref{alg:Algm}\ with $m=1$, assume that there exists $\beta>0$
such that $\Bdivx{h}(f,f') \geq \beta \Bdivx{\loss_y}(f,f')$ for any $f$ and $f'$, and 
assume that $\alpha \in (0,\beta]$.  
Assume also that $Q_t(\theta)$ defined in Algorithm \ref{alg:Algm}\ is $1/\eta$ smooth in $\theta$:
\[
\|\nabla Q_t(\theta) - \nabla Q_t(\theta')\| \leq (1/\eta) \|\theta - \theta'\| .
\]

Assume that $\theta_{t+1}$ is an improvement of $\theta_t$ with respect to minimizing $Q_t$ so that 
\[
 Q_t(\theta_{t+1}) \leq Q_t(\tilde{\theta}) , 
\]
where
\begin{align*}
\tilde{\theta} = \theta_t - \eta \nabla Q_t(\theta_t) . 
\end{align*}
Then we have 
\[
  \reglossx{\alpha}(\theta_{t+1}) \le \reglossx{\alpha}(\theta_t) - \frac{\alpha \eta}{2}
  \|\nabla \reglossx{\alpha}(\theta_{t})\|^2 .
\]
}
\begin{proof}

We first define $\tilde{Q}_t(\theta)$ as follows:
\[
\tilde{Q}_t(\theta) := \meanOverS{\Bdivx{h}(\netprm,\netprmt) + \alpha \nabla
  \loss_y(\netprmt)^\top \netprm } + R(\theta) .
\]
We can check that $Q_t(\theta) - \tilde{Q}_t(\theta)$ is independent of
$\theta$. Therefore optimizing $\theta$ with respect to $Q_t(\theta)$ is the
same as optimizing $\theta$ with respect to $\tilde{Q}_t(\theta)$, and
$\nabla Q_t(\theta)=\nabla \tilde{Q}_t(\theta)$.

The smoothness assumption implies that
\[
\tilde{Q}_t(\theta- \Delta \theta) \leq \tilde{Q}_t(\theta) - \nabla {Q}_t(\theta)^\top \Delta
\theta + \frac{1}{2\eta } \|\Delta \theta\|^2 .
\]
Therefore
\begin{align*}
\tilde{Q}_t(\theta_{t+1}) \leq& \tilde{Q}_t(\tilde{\theta})
= \tilde{Q}_t(\theta_t - \eta \nabla Q_t(\theta_t)) \\
\leq & \tilde{Q}_t(\theta_t) - \eta \|\nabla Q_t(\theta_t)\|^2 +
       \frac{1}{2\eta} \| \eta \nabla Q_t(\theta_t)\|^2 \\
=& \tilde{Q}_t(\theta_t) - \frac{\eta}{2} \|\nabla Q_t(\theta_t)\|^2 .
\end{align*}
Note also that
\begin{align*}
\tilde{Q}_t(\theta_{t+1}) - \tilde{Q}_t(\theta_t) \geq&
\meanOverS{\beta \Bdivx \loss_y(f_{\theta_{t+1}},f_{\theta_t}) + \alpha \nabla \loss_y(f_{\theta_t})^\top
                                 (f_{\theta_{t+1}} -f_{\theta_t}) }
+ [R(\theta_{t+1}) - R(\theta_t)]\\
=&
\meanOverS{ 
(\beta-\alpha) \Bdivx \loss_y(f_{\theta_{t+1}},f_{\theta_t}) + \alpha
   \loss_y(f_{\theta_{t+1}}) - \alpha \loss_y(f_{\theta_t})}
+ [R(\theta_{t+1}) - R(\theta_t)]\\
\geq&
\meanOverS{ 
\alpha   \loss_y(f_{\theta_{t+1}}) - \alpha \loss_y(f_{\theta_t})}
+ [R(\theta_{t+1}) - R(\theta_t)]\\
=& \alpha \reglossx{\alpha}(\theta_{t+1}) - \alpha \reglossx{\alpha}(\theta_{t}) .
\end{align*}
The second inequality is due to the non-negativity of the Bregman divergence.  

By combining the two inequalities, we obtain
\[
\alpha \reglossx{\alpha}(\theta_{t+1}) \leq \alpha \reglossx{\alpha}(\theta_{t}) -
\frac{\eta}{2} \|\nabla Q_t(\theta_t)\|^2 .
\]
Now, observe that $\nabla Q_t(\theta_t) = \nabla
\tilde{Q}_t(\theta_t) =\alpha \nabla
\reglossx{\alpha}(\theta_t)$, and we obtain the desired bound.
\end{proof}


\section{Empirical study} 
\label{sec:empirical}
While the proposed framework is general, our empirical study places a
major focus  on \oursL\ (Algorithm \ref{alg:AlgL}) with the cross-entropy loss, due to its connection to distillation 
(Proposition \ref{pr:distill}).  
In particular, we set up our implementation so that one instance of \oursL\ 
coincides with 
self-distillation 
to provide empirical insight into it from a functional gradient viewpoint.  

First, with the goal of understanding the empirical behavior of the algorithm,  
we examine obtained models in reference to our theoretical findings.  
We use relatively small neural networks for this purpose.  
Next, we study the case of larger networks with consideration of practicality.  

\subsection{Implementation}
\label{sec:implement}

To implement the algorithms presented above, methods of 
parameter initialization and optimization need to be considered. 
To observe the basic behavior, our strategy in this work is to keep it as simple as possible.  

\tightpara{Initial parameter $\theta_0$}
As the functions of interest are nonconvex, the outcome depends on the initial parameter $\theta_0$.  
The most natural (and simplest) choice is random parameters.  
This option is called `ini:random' below.  
We also considered two more options. %
One is to start from a {\em base model} obtained by regular training, 
called `ini:base'.  
This option enables study of self-distillation.  
The other is to start from a shrunk version of the base model, 
and details of this option will be provided later.  %

\tightpara{Parameter update}
To update parameter $\theta$ by descending the stochastic gradient, 
standard techniques can be used such as momentum, Rmsprop \cite{rmsprop}, 
Adam \cite{Adam15}, and so forth.  
As is the case for regular training, learning rate scheduling is beneficial.  
Among many possibilities, we chose to repeatedly use for each $t$, 
the same method that works well for regular training.  
For example, a standard method for CIFAR10 is 
to use momentum and decay the learning rate only a few times, 
and therefore we use this scheme for each stage on CIFAR10.  
That is, the learning rate is reset to the initial rate for each $t$; 
however note that $\theta$ is not reset.  
Although this is perhaps not the best strategy in terms of computational cost, 
its advantage is that at the end of each stage, we obtain ``clean'' intermediate models 
with $\theta_t$ that were optimized for intermediate goals. 
(If instead, we used one decay schedule from the beginning to the end, 
the convergence theorem still holds, but 
$\theta_t$ would be noisy when the learning rate is still high.) 
This strategy enables to study how a model changes as the guide function gradually goes ahead, 
and also relates the method to self-distillation. 


Since $\theta$ is not reset when the learning rate is reset, this schedule can be regarded as 
a simplified fixed-schedule version of {\em SGD with warm restarts (SGDR)} \cite{cosine17}. 
(SGDR instead does sophisticated scheduling with cosine-shape decay and variable epochs.) 
For comparison, we test the same schedule with the standard optimization objective 
(`base-loop'; Algorithm \ref{alg:baseLoop}). 

\tightpara{Enabling study of self-distillation}
We study classification tasks with the standard cross-entropy loss, which satisfies the condition 
of Proposition \ref{pr:distill}.  
Combined with the choice of learning rate scheduling above, 
\oursL\ with the ini:base\ option (which initializes $\theta_0$ with a trained model) 
essentially becomes self-distillation. 
Thus, one aspect of our experiments is 
to study self-distillation from the viewpoint of functional gradient learning.   

\begin{algorithm}[t]
\begin{small}
   \caption{base-loop\ (simplified SGDR): 
         {\bfseries Input:} $\theta_0$, training set $S$.  
         Meta-parameter: $ T $. 
         {\bfseries Output:} $\theta_ T $. 
   } 
   \label{alg:baseLoop}
\begin{algorithmic}
   \FOR{$t=0$ {\bfseries to} $ T -1$}  
      \STATE $\theta_{t+1}  \leftarrow  \arg\min_\theta \left[ 
               \meanOverSp{ \loss_y(\net{\theta}{x}) } + R(\theta)
             \right]$ \\
             $~~~~~~~~~~~~~$ where $\theta$ is initialized by $\theta_t$.  
   \ENDFOR
\end{algorithmic}
\end{small}
\end{algorithm}

\subsection{Experimental setup}
\label{sec:setup}
\begin{table}[t]
\begin{center} \begin{small} \begin{tabular}{|c|r|r|r|r|}

\hline
           & \#class & \multicolumn{1}{|c|}{train} & \multicolumn{1}{|c|}{dev.} & \multicolumn{1}{|c|}{test} \\
\hline
  CIFAR10  &   10 &   49000 &  1000 & 10000 \\
  CIFAR100 &  100 &   49000 &  1000 & 10000 \\
  SVHN     &   10 &  599388 &  5000 & 26032 \\
  ImageNet & 1000 & 1271167 & 10000 & 50000 \\
\hline
\end{tabular} \end{small} \end{center}
\vskip -0.15in
\caption{\label{tab:data-stat}
  Data. 
For each dataset, we randomly split the official training set into 
a training set and a development set to use
the development set for meta-parameter tuning.  
For ImageNet, following custom, %
we used the official validation set as our `test' set.   
}
\vskip -0.1in
\end{table}
Table \ref{tab:data-stat} summarizes the data we used.  
As for network architectures, we mainly used ResNet 
\cite{Resnet16a,Resnet16b} and wide ResNet (WRN) \cite{Wresnet16}. 
Following the original work, the regularization term $R(\theta)$ was 
set to be $R(\theta)= \frac{\lambda}{2}\| \theta \|^2$ 
where $\lambda$ is the weight decay. 
We fixed mini-batch size to 128 and used the same learning rate decay schedule 
for all but ImageNet.  
Due to the page limit, details are described in the supplementary material.  
However, 
note that the schedule we used for all but ImageNet
is 3--4 times longer than those used in the original ResNet or WRN study
for CIFAR datasets.  
This is because we used the ``train longer'' strategy \cite{HHS17}, and 
accordingly, the base model performance visibly improved from the original work.  
This, in fact, made it harder to obtain large performance gains 
over the base models 
(not only for GULF\ but also for all other tested methods) 
as the bar was set higher.  We feel that this is more realistic 
testing than using the original shorter schedule.  

We applied the standard mean/std normalization to images 
and used the standard image augmentation.  
In particular, for ImageNet, we used the same data augmentation scheme 
as used for training the pre-trained models provided as part of TorchVision, 
since we used these models as our base model.  

The default value of $\alpha$ is 0.3. 

\begin{figure}[h]
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-with-legend}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-with-legend}
  \oursL
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-a0-tloss-with-legend}}
\vskip -0.1in
\caption{ \label{fig:c100k1-a0-tloss-with-legend}
  Regular training
}
\end{center}
\end{subfigure}%
\vskip -0.3in
\caption{ \label{fig:c100k1-tloss}
Test loss in relation to training loss. 
The arrows indicate the direction of time flow. 
CIFAR100.  ResNet-28.
} 
\vskip -0.1in
\end{figure}

\begin{figure*}
\centering
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-a1}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-a1}
  base-loop\ ($\alpha$=1)
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-a09}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-a09}
  $\alpha$=0.9
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-a03}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-a03}
  $\alpha$=0.3
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-a01}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-a01}
  $\alpha$=0.1
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.2\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-tloss-a001}}
\vskip -0.1in
\caption{ \label{fig:c100k1-tloss-a001}
  $\alpha$=0.01
}
\end{center}
\end{subfigure}%
\vskip -0.1in
\vskip -0.25in
\caption{ \label{fig:c100k1-tloss-alpha}
Test loss %
of ini:base (`$\triangle$') and ini:random(`$\circ$').  
with five values of $\alpha$ (becoming smaller from left to right), 
in relation to training loss. \oursL.  $ T $=25. 
CIFAR100.  ResNet-28. 
As $\alpha$ becomes smaller, 
the (potential) meeting point shifts 
further away from the base model. 
The left-most figure shows base-loop, which is equivalent to 
$\alpha$=1. 
}
\vskip -0.1in
\end{figure*}



\subsection{Smooth path} 

We start with examining training of a relatively small network ResNet-28 (0.4M parameters) 
on CIFAR100.  
In this setting, optimization is fast, and so a relatively large $ T $ 
(the number of stages) is feasible.  

We performed \oursL\ training with $ T $=25 starting from random parameters (ini:random)
as well as starting from a base model obtained by regular training (ini:base). 
Figure \ref{fig:c100k1-tloss-with-legend} plots test loss of these two runs 
in relation to training loss. 
Each point represents a model $\net{\theta_t}{x}$ 
at time $t=1,3,5,\cdots,25$, and the arrows
indicate the direction of time flow.  
We observe that training proceeds on a {\em smooth path}.  
ini:random ($\circ$), which starts from random parameters ($\square$), 
reduces both training loss and test loss.  
ini:base ($\triangle$) starts from the base model ($\times$) and 
increases training loss, but reduces test loss.  
ini:random\ and ini:base\ meet and complete one {\em smooth path} from a random state ($\square$)
to the base model ($\times$).  
ini:random\ goes forward on this path while ini:base\ goes backward, 
and importantly, the path goes through the region where test loss is lower than 
that of the base model.  
The test error plotted against training loss also forms a U-shape path.  
Similar U-shape curves %
were observed across datasets and network architectures. 
The supplementary material shows 
a test error curve and 
a few more examples of test loss curves including a case of DenseNet \cite{densenet17}. 
 
In the middle of this path, a number of models with good generalization performance lie.  
One might wonder if regular training also forms such a path.   
Figure \ref{fig:c100k1-a0-tloss-with-legend} shows that this is not the case. 
This figure plots the loss of intermediate models 
in the course of regular training so that 
the $i$-th point represents a model after 20K$\times i$ steps of 
mini-batch SGD with the learning rate being reduced twice.   
The path of regular training from random initialization ($\square$) to the final model ($\times$) 
is rather bumpy and the test loss generally stays as high as the final outcome.   
The bumpiness is due to the fact that the learning rate is relatively high at the beginning of training.  
Comparing Figures \ref{fig:c100k1-tloss-with-legend} and \ref{fig:c100k1-a0-tloss-with-legend}, 
GULF\ training clearly takes a very different path from regular training.  

\subsection{In relation to the theory} 

\tightpara{Going forward, going backward} 
It might look puzzling why ini:base\ goes {\em backward} in the direction of {\em increasing} the training loss.  
Theorem \ref{thm:L-conv} suggests that this is the effect of the regularization term $R(\theta)$, 
in this case $R(\theta)= \frac{\lambda}{2}\|\theta\|^2$ with weight decay $\lambda$. 
The theory indicates that for $\alpha \in (0,1]$, the $\alpha$-regularized loss 
\[
  \reglossx{\ssa}(\prm) = \meanOverS{ \loss_y(\net{\theta}{x}) } + R(\theta)/\alpha
\]  
goes down and eventually converges as \oursL\ proceeds.  By contrast, 
The base model is a result of minimizing 
\[
  \meanOverS{ \loss_y(\net{\theta}{x}) } + R(\theta). %
\]
As we always set $\alpha<1$ (0.3 in this case), i.e., $1/\alpha>1$, 
\oursL\ prefers smaller parameters than the base model does.  
Consequently, 
when \oursL\ (with small $\alpha$) starts from the base model (which has low training loss and high $R(\theta)$ ), 
\oursL\ is likely to reduce $R(\theta)/\alpha$ at the expense of increasing loss (going backward).  
When \oursL\ starts from random parameters, whose training loss is high, 
\oursL\ is likely to reduce loss (going forward) 
at the expense of increasing $R(\theta)/\alpha$. 

\tightpara{Effects of changing $\alpha$} 
With \oursL, the guide function $ \netfunc^* $ satisfies
\begin{align*}
   \netfunc^*  &\approx \netprmt - \alpha (\hessian{\loss_y(\netprmt)})^{-1} \nabla \loss_y(\netprmt), 
\end{align*}
thus, $\alpha$ serves as a step-size of functional gradient descent for {\em reducing loss}.  
The effects of changing $\alpha$ are shown in Figure \ref{fig:c100k1-tloss-alpha} 
with $ T $ fixed to 25.  
The left-most graph is base-loop, which is equivalent to \oursL\ with $\alpha$=1 in 
this implementation. 
There are three things to note.  
First, with a very small step-size $\alpha$=0.01 (the right most), ini:random\ cannot reach far from the random state 
for $ T $=25. 
This is a straightforward effect of a small step size.   
Second, 
as step-size $\alpha$ becomes smaller (from left to right), 
the (potential) meeting/convergence point shifts further away from the base model; 
the convergence point of $\|\theta_t\|^2$ also shifts away from the base model 
and decreases (supplementary material). 
This is the effect of larger $R(\theta)/\alpha$ for smaller $\alpha$.  
Finally, with a large step-size (0.9 and 1), the curve flattens and 
it no longer goes through the high-performance regions {\em slowly} or {\em smoothly}, 
and the benefit diminishes/vanishes. 

\tightpara{$\alpha$-regularized loss $\reglossx{\ssa}(\prm)$} 
Figure \ref{fig:c100k1-ajl2l-alpha} confirms that, as suggested by the theory, 
 $\reglossx{\ssa}(\prm)$ goes down and almost converges as training proceeds.  
This fact motivates examining standard models trained with this $\reglossx{\ssa}(\prm)$ objective, 
which we call base-$\lambda/\ssa$\ models.  
We found that base-$\lambda/\ssa$\ models do not perform as well as \oursL\ at all.  In particular, 
with a very small $\alpha$=0.01, which 100 times tightens regularization, 
test error of base-$\lambda/\ssa$\ drastically degrades due to underfitting; 
in contrast, ini:base\ with $\alpha$=0.01 performs well. 
Moreover,  
base-$\lambda/\ssa$\ models are very different from \oursL\ models with corresponding $\alpha$
even with a moderate $\alpha$.  
For example, 
Figure \ref{fig:c100k1-l2-baselam-a03} plots the parameter size $\|\theta_t\|^2$ 
in relation to training loss for $\alpha$=0.3. base-$\lambda/\ssa$\ is clearly far away from 
where ini:base\ and ini:random\ converge to.  




\begin{figure}[h]
\vskip -0.1in
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-ajl2l-alpha}}
\vskip -0.1in
\caption{ \label{fig:c100k1-ajl2l-alpha}
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k1-l2-baselam-a03-with-legend}}
\vskip -0.1in
\caption{ \label{fig:c100k1-l2-baselam-a03}
}
\end{center}
\end{subfigure}%
\vskip -0.3in
\caption{ \label{fig:c100k1-ajl2l-baselam-alpha}
(a) $\alpha$-regularized loss $\reglossx{\ssa}(\prm)$ in relation to time $t$. \oursL\ ini:base. 
(b) $\|\theta_t\|^2$ and training loss of base-$\lambda/\ssa$\ in comparison with \oursL.  $\alpha$=0.3. 
CIFAR100.  ResNet-28. 
}
\end{figure}
\vskip -0.1in


\tightpara{Benefit of guiding} 
This fact illustrates the merit of guided learning (including self-distillation). 
GULF\ (indirectly and locally) minimizes the $\alpha$-regularized loss $\reglossx{\ssa}(\prm)$, but it does this 
against the restraining force of {\em pulling the model back} to the current model.  This serves as 
a form of regularization.  Without such a force, training for, say, $\reglossx{0.01}$ 
would make a big jump to rapidly reduce the parameter size and end up with 
a radical solution that suffers severely from underfitting.  This is what happens with base-$\lambda/\ssa$. 
By contrast, guided learning finds a more moderate solution with good generalization performance, 
and this is the benefit of extra regularization (in the form of pulling back) provided through the 
guide function.  
The regularization effect of distillation has been mentioned \cite{distill14}, 
and our framework formalizes the notion through the functional gradient learning viewpoint. 

\subsection{With smaller networks}

\newcommand{\imp}{\em} %
\begin{table}[b]
\vskip -0.2in
\begin{center} \begin{small} \begin{tabular}{c|c|c|c|c|c|c|}
\cline{2-7}
  &\multicolumn{2}{|c|}{}& C10 & C100 & \multicolumn{2}{|c|}{SVHN} \\
\cline{2-7}
1 & \multirow{4}{*}{baselines} 
  & base model  & 6.42 & 30.90 & 1.86 & 1.64 \\ 
\cline{3-7}  
2 && base-$\lambda/\ssa$     & 6.60 & 30.24 & 1.78 & 1.67     \\ %
3 && base-loop    & 6.20 & 30.09 & 1.93 &{\bf 1.53}\\
4 && label smooth & 6.66 & 30.52 & 1.71 & 1.60     \\
\cline{2-7}
5 & \multirow{2}{*}{\oursL}   
   & ini:random &     5.91 &{\bf 28.83} &     1.71 &{\bf 1.53} \\
6 && ini:base   &{\bf 5.75}&     29.12  &{\bf 1.65}&     1.56  \\
\cline{2-7}
\end{tabular} \end{small} \end{center}
\vskip -0.125in
\caption{\label{tab:c10-c100-svhn}
  Test error (\%).  Median of 3 runs.
  Resnet-28 (0.4M parameters) for CIFAR10/100, 
  and WRN-16-4 (2.7M parameters) for SVHN.  
  Two numbers for SVHN are without and with dropout. 
  base-$\lambda/\ssa$: weight decay $\lambda/\alpha$.  
  base-loop: Algorithm \ref{alg:baseLoop}.  
}
\vskip -0.1in 
\end{table}


Now we review the test error results of using relatively small networks 
in Table \ref{tab:c10-c100-svhn}.  
$ T $ for \oursL\ and base-loop\ was fixed to 25 on CIFAR10/100 
and 15 on SVHN.  
Step-size $\alpha$ was fixed to 0.3 for ini:random\ and 
chosen from $\{$ 0.01,0.03 $\}$ for ini:base.  
\oursL\ is consistently better than the base model (Row 1)
and generally better than the three baseline methods (Row 2--4).  
The base-$\lambda/\ssa$\ results (Row 2) were obtained by $\alpha$=0.3, and they are generally 
not much different from the base model.  
base-loop\ (Row 3) generally makes small improvement over the base model, 
but it generally falls short of \oursL.  
A common technique, label smoothing (Row 4) \cite{inception16}, 
`softens' labels by taking a small amount of 
probability from the correct class and distributing it equally to the incorrect classes.  
It 
generally worked well, but the improvements were small. 
That is, 
the three baseline methods produced performance gains to some extent, but 
their gains are relatively small, and they are not as consistent as \oursL\ across datasets.  

\tightpara{ini:random} 

In these experiments, 
ini:random\ performed as well as ini:base.  This fact cannot be explained from 
the knowledge-transfer viewpoint of distillation, but it can be explained from our 
functional gradient learning viewpoint, as in the previous section.  

\begin{figure}[t]
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c10k10-tloss-with-fcs}}
\vskip -0.1in
\caption{ \label{fig:c10k10-tloss-with-fcs}
  CIFAR10
}
\end{center}
\end{subfigure}%
\begin{subfigure}[b]{0.5\linewidth}
\begin{center}
\centerline{\includegraphics[width=1\linewidth]{aug-c100k10-tloss-with-fcs}}
\vskip -0.1in
\caption{ \label{fig:c100k10-tloss-with-fcs}
  CIFAR100
}
\end{center}
\end{subfigure}%
\vskip -0.25in
\caption{ \label{fig:c10c100k10-tloss-with-fcs}
  Test loss in relation to training loss. 
  WRN-28-10 on CIFAR10 and CIFAR100. \oursL. 
  ini:base/2 (`$\diamond$') fills the gap between ini:random\ (`$\circ$') and ini:base\ (`$\triangle$'). 
}
\end{figure}

\begin{table}[t]
\begin{center} \begin{small} \begin{tabular}{c|c|c|c|}
\cline{2-4}
&     & CIFAR10 & CIFAR100 \\
\cline{2-4}
1&base model& 3.82 & 18.55 \\ 
\cline{3-4}  
2&base-$\lambda/\ssa$    & 3.70 & 27.89 \\
3&base-loop  & 3.70 & 18.91 \\  %
4&lab smooth & 4.13 & 19.44 \\  %
\cline{2-4}
5& \ourssq  &{\bf 3.46}& 18.14 \\ %
6& \oursL   & 3.63 &{\bf 17.95}\\ %
\cline{2-4}
\end{tabular} \end{small} \end{center}
\vskip -0.15in
\caption{\label{tab:c10-c100-k10}
Test error (\%) results on CIFAR10 and CIFAR100.  
WRN-28-10 (36.5M parameters) without dropout.  
Median of 3 runs. 
}
\vskip -0.1in
\end{table}

\subsection{With larger networks} 
The neural networks and the size of images (32$\times$32) used above are relatively small.  
We now consider computationally more expensive cases. %

\tightpara{Parameter shrinking}
ini:random, the most natural option from the functional gradient learning viewpoint, 
unfortunately, turned out to be too costly in this large-network situation.  
Moreover, in this setting, 
it is useful to have an option of starting 
somewhere between the two end points (`random' and `base') 
since that is where good models tend to lie according to our study with small networks.  
Therefore, we experimented with `rewinding' a base model 
by shrinking its weights and bias of the last fully-connected linear layer 
by dividing them with $\fcsmeta > 1$ (a meta-parameter).  
We use these partially-shrunk parameters as the initial parameter $\theta_0$ for 
GULF.  Since doing so shrinks the model output $\net{\theta_0}{x}$ 
by the factor of $\fcsmeta$, this is closely related to {\em temperature scaling}, 
for {\em distillation} \cite{distill14} and post-training calibration \cite{GPSW17}.  
Parameter shrinking is, however, simpler than temperature scaling of distillation, 
which scales logits of both models, 
and fits well in our framework. 

Figure \ref{fig:c10c100k10-tloss-with-fcs} shows 
training loss (the $x$-axis) and test loss (the $y$-axis) 
obtained when parameter shrinking is applied to WRN-28-10 
on CIFAR10 and CIFAR100.  
By shrinking with $\fcsmeta$=2, the loss values of the base model change 
from `base' ($\times$)
to `base/2' ($+$). 
The location of base/2 is roughly the midpoint of two end points 
`base' and `random'. 
ini:base/2 ($\diamond$), which starts from the shrunk model, explores the space
neither ini:random\ nor ini:base\ can reach in a few stages. 

\tightpara{Larger ResNets on CIFAR10 and CIFAR100}

Table \ref{tab:c10-c100-k10} shows test error of ini:base/2 
using WRN-28-10 on CIFAR10/100. 
$ T $ was fixed to 1.  
Compared with the base model, 
both \ours1 and 2 consistently improved performance, 
while the baseline methods mostly failed to make improvements.  
\ours1 and 2 produced similar performances. %
This is the best WRN non-ensemble results on CIFAR10/100 
among the self-distillation related studies that we are aware of.  



\tightpara{ImageNet}

To test further scale-up, 
we experimented with 
ResNet-50 (25.6M parameters) and WRN50-2 (68.9M parameters)
on the ILSVRC-2012 ImageNet dataset. 
As ImageNet training is resource-consuming, we only tested 
selected configurations, which are \oursL\ with ini:base\ and ini:base/2 options. 
In these experiments, $\alpha$ was set to 0.5, but partial results suggested that 0.3  
works well too.  
We used models pre-trained on ImageNet provided as part of TorchVision\footnote{
  https://pytorch.org/docs/stable/torchvision/models.html
}
as the base models.  
Table \ref{tab:imgnet} shows that 
\oursL\ consistently improves error rates over the base model.  
The best-performing ini:base/2 achieved lower error rates than 
a twice deeper counterpart of each network, 
ResNet-101 for ResNet-50 and WRN-101-2 for WRN-50-2
trained in a standard way (Rows 11--12).  
Thus, we confirmed that \oursL\ scales up and brings performance gains on ImageNet.  
To our knowledge, this is one of the largest-scale 
ImageNet experiments among the self-distillation related studies.  

\begin{table}[h]
\begin{center} \begin{small} \begin{tabular}{c|c|c|cc|cc|}
\cline{2-7} %
&\multicolumn{2}{|c|}{methods}
                   & \multicolumn{2}{|c|}{Resnet-50}
                   & \multicolumn{2}{|c|}{WRN50-2} \\                  
\cline{2-7} \noalign{\vskip\doublerulesep \vskip-\arrayrulewidth} \cline{2-7}
1&\multicolumn{2}{|c|}{base model} & 23.87 & 7.14 & 21.53 & 5.91 \\ %
\cline{2-7}
2&\multirow{3}{*}{base-loop}
           & $t$=1 & 23.73 & 6.95 & 21.99 & 6.11 \\
3&         & $t$=2 & 23.50 & 6.93 & \multicolumn{2}{|c|}{--} \\
4&         & $t$=3 & 23.36 & 6.78 & \multicolumn{2}{|c|}{--} \\
\cline{2-7} \noalign{\vskip\doublerulesep \vskip-\arrayrulewidth} \cline{2-7}
5&\multirow{3}{*}{ini:base}
           & $t$=1 & 22.79 & 6.43 & 21.17 & 5.65 \\ %
6&         & $t$=2 & 22.49 & 6.27 & \multicolumn{2}{|c|}{--} \\
7&         & $t$=3 & 22.31 & 6.28 &       &      \\    
\cline{2-7}
\cline{2-7}
8&\multirow{3}{*}{ini:base/2}          
           & $t$=1 &     22.50  &     6.25 &{\bf 20.69}&{\bf 5.35}\\
9&         & $t$=2 &     22.31  &     6.18 & \multicolumn{2}{|c|}{--}  \\
10&         & $t$=3 &{\bf 22.08} &{\bf 6.10}& \multicolumn{2}{|c|}{--}  \\
\cline{2-7} \noalign{\vskip\doublerulesep \vskip-\arrayrulewidth} \cline{2-7}
11&\multicolumn{2}{|c|}{Resnet-101$\dagger$} & 22.63 & 6.44 & \multicolumn{2}{|c|}{--} \\
12&\multicolumn{2}{|c|}{WRN-101-2$\dagger$}  &  \multicolumn{2}{|c|}{--} & 21.16 & 5.72 \\
\cline{2-7}
\end{tabular} \end{small} \end{center}
\vskip -0.15in
\caption{\label{tab:imgnet}
ImageNet 224$\times$224 single-crop results on the validation set. 
\oursL. 
top-1 and top-5 errors (\%). \\
$\dagger$
The ResNet-101 and WRN-101-2 performances are from the description of 
  the pre-trained torchvision models.  
}
\vskip -0.1in
\end{table}


\tightpara{Additional experiments on text}
Finally, the experiments in this section used image data.  Additional experiments 
using text data are presented in the supplementary material.  


\section{Discussion}

\tightpara{Guided exploration of landscape}
GULF\ is an informed/guided exploration of the loss landscape, 
where the guidance is successively given as interim goals set in the neighborhood of the model at the time, 
and such guidance is provided by gradient descent in a function space.  
Another view of this process is an accumulation of successive greedy optimization.  
Instead of searching the entire space for the ultimate goal of loss minimization at once, 
guided learning proceeds with repeating a local search, which limits the space to be searched
and leads to better generalization.  
Its benefit is analogous to that of $\varepsilon$-boosting. 


\tightpara{\ourssq}
\oursL\ uses the second-order information of loss in the functional gradient step 
for generating the guide function, and \ourssq\ does not.  
\oursL's update rule is equivalent to that of distillation, and 
\ourssq's is not. \ourssq\ also differs from the logit least square fitting version of distillation.  
In our experiments (though limited due to our focus on self-distillation study), 
\ourssq\ performed as well as \oursL. %
If this is a general trend, this indicates that 
inclusion of the second-order information is not particularly helpful.  
If so, this could be because the second-order information is useful for accelerating optimization, 
but we would like to proceed slowly to obtain better generalization performance.  
This motivates further investigation of \ourssq\ as well as other instantiations of the framework.  

\tightpara{Computational cost}
From a practical viewpoint, a shortcoming of the particular setup tested here 
(but not the general framework of GULF) 
is computational cost.  
Since we used the same learning rate scheduling as regular training in each stage, 
GULF\ training with $ T $ stages took more than $ T $ times longer than 
regular training.  
It is conceivable that training in each stage can be shortened without hurting 
performance since optimization should be easier as a results of aiming at a nearby goal. 
Schemes that decay the learning rate throughout the training without restarts 
or 
hybrid approaches
might also be beneficial 
for reducing computation.  Note that Theorem \ref{thm:L-conv} does not require each 
stage to be performed to the optimum.  
On the other hand, testing (i.e., making predictions) of the models trained with GULF\ 
only requires the same cost as regular models.  
As shown in the ImageNet experiments, a model trained with GULF\ 
could perform better than a much larger (and so slow-to-predict) model; 
in that case, GULF\ can save the overall computational cost 
since the cost for making predictions can be significant 
for practical purposes.  



\tightpara{Relation to other methods} 
The proposed method seeks to improve generalization performances in a principled way 
that limits the searched parameter space.  The relation to existing methods for similar purposes 
is at least two-fold.  First, we view that this work gives theoretical insight into related methods 
such as self-distillation and label smoothing, which we hope can be used to improve them.  
Second, methods derived from this framework can be used {\em with} existing techniques 
that are based on different principles (e.g., weight decay and dropout) for {\em further} improvements.  

\tightpara{Distillation} 
Due to the connection discussed above, 
our theoretical and empirical analyses of \oursL\ provide a new functional gradient view of distillation.  
Here we discuss a few self-distillation studies from this new viewpoint. 
\cite{Ban18} showed that iterative self-distillation improves performance over the base model. 
They set $\alpha$ to 0 (in our terminology) and 
reported that there were no performance gains on CIFAR10. 
According to our theory, 
when $\alpha$ goes to 0, the quantity reduced throughout the process is 
not the $\alpha$-regularized loss but merely $R(\theta)$.  
Such an extreme setting might be risky.  
In {\em deep mutual learning} \cite{deepMutual18}, 
multiple models are simultaneously trained by reducing loss and aligning 
each other's model output.  They were surprised by the fact that 
`no prior powerful teacher' was necessary. 
This fact can be explained by our functional gradient view by 
relating their approach to our ini:random.  
Finally, the regularization effect of distillation has been noticed \cite{distill14}. 
Our framework formalized the notion through the functional gradient learning viewpoint. 





\subsection{Additional experiments on text data}
We tested GULF\ on sentiment classification
to predict whether reviews are positive or negative, 
using the polarized Yelp dataset (\#train: 560K, \#test: 38K) \cite{ZZL15}.  
The best-performing models on this task are transformers pre-trained with language modeling on 
large and general text data such as Bert \cite{bert19} and XLnet \cite{xlnet19}.  
However, these models are generally large and time-consuming to train
using a GPU (i.e., without TPUs used in the original work).  
Therefore, instead, we used the deep pyramid convolutional neural network (DPCNN) 
\cite{dpcnn17} as our base model. 
In these experiments, we used \oursL.  

Table \ref{tab:text-results} shows the test error results in five settings.  
The last three use relatively small training sets of 45K data points 
and validation sets of 5K data points, randomly chosen from the original training set, 
while the first two use the entire training set (560K data points) except for 
5K data points held out for validation (meta-parameter tuning).  
DPCNNs optionally take additional features produced by embeddings of text regions that are trained with unlabeled data, 
similar to language modeling.  
Cases 1 and 3 exploited this option, training embeddings using the entire training set as unlabeled data; 
\ref{apx:text-details} below provides the details.  
As in the image experiments, we used the cross entropy loss with softmax except for 
Case 5, where the quadratic hinge loss $L_y(f) = \max(0,1-yf)^2$ for $y \in \{-1,1\}$ was used.  
This serves as an example of extending self-distillation (formulated specifically with the cross-entropy loss)
to general loss functions. 

In all the five settings, GULF\ achieves better test errors than the baseline methods, 
which shows the effectiveness of our approach in these settings.  
On this task, dropout turned out to be not very effective, 
which is, however, a reminder that the effectiveness of regularization methods can be data-dependent in general. 

 
\begin{table}[h]
\vskip -0.1in
\begin{center} \begin{small} \begin{tabular}{|c|l|c|c|c|c|}
  \multicolumn{2}{r}{Case\#}
 & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{2} & \multicolumn{1}{c}{} \\
\hline
\multicolumn{2}{|c|}{}
  &\multicolumn{2}{|c|}{LM-like prep?}                        & Runtime &Text for \\ %
\cline{3-4}
\multicolumn{2}{|c|}{}
  &\multicolumn{1}{|c|}{$~\;$Yes$~\;$}&\multicolumn{1}{|c|}{$~\;$No$~\;$} & (sec/K)   &prep (GB)\\ %
\hline
(J \& Z, 2017)&DPCNN       &{\em 2.64}&{\em 3.30}&$~$0.1 &0.4\\%& 10M \\
\hline
\multirow{2}{*}{This work}
&Table \ref{tab:text-results} best& 2.34 & 2.70 &$~$0.1 &0.4 \\%&10M \\ %
&Ensemble                 &{\bf 2.18}&{\bf 2.46}&$~$0.9 &0.4 \\%& 10M \\ %
\hline
\multirow{2}{*}{\cite{bert19}}
                &Bert base &     2.25 &    6.19 &$~$5.7        &13 \\%& 110M \\ %
                &Bert large&{\em 1.89}&   --    & 17.9         &13  \\%& 335M \\ %
\hline
\multirow{2}{*}{\cite{xlnet19}}
                &XLnet base&     1.92 &    4.51 & 17.2         &13  \\%& 117M \\ %
                &XLnet large&{\em 1.55}&  --    & 40.5         &126 \\%& 361M \\ %
\hline
\end{tabular} \end{small} \end{center}
\vskip -0.1in
\caption{\label{tab:text-results2}
  GULF\ ensemble results on Yelp in comparison with previous models.
  Test error (\%) with or without embedding learning (DPCNN) or 
  language modeling-based pre-training (Bert and XLnet), respectively, 
  corresponding to Cases 1 \& 2 of Table \ref{tab:text-results}.  
  Runtime: real time in seconds for labeling 1K instances using a single GPU with 11GB device memory, 
  measured in the setting of Case 1; the average of 3 runs.  
  The last column shows 
  amounts of text data in giga bytes used for pre-training or embedding learning in Case 1.  \\
  The test errors in {\em italics} were copied from the respective publications 
  except that the Bert-large test error is from \cite{XDHLL19}; 
  other test errors and runtime were obtained by our experiments.  
  Our ensemble test error results are in bold.  
}
\end{table}

It is known that performance can be improved by making an ensemble of models from different stages 
of self-distillation, e.g., \cite{Ban18}.  
In Table \ref{tab:text-results2}, we report ensemble performances of 
DPCNNs trained with GULF, in comparison with the previous best models.  
Test errors with and without embedding learning 
(or language modeling-based pre-training for Bert and XLnet) are shown, 
corresponding to Cases 1 and 2 in Table \ref{tab:text-results}.  
The ensemble results were obtained by adding after applying softmax 
the output values of 20 DPCNNs (or 10 in Case 2) of last 5 stages 
of GULF\ training with different training options; details are provided in \ref{apx:text-details}. 

With embedding learning, the ensemble of DPCNNs trained with GULF\ 
achieved test error 2.18\%, which slightly beats 2.25\% of pre-trained Bert-base, 
while testing (i.e., making predictions) of this ensemble is more than 6 times faster than Bert-base, 
as shown in the `Runtime' column.  
(Note, however, that runtime depends on implementation and hardware/software configurations.)
That is, using GULF, we were able to obtain a classifier that is as accurate as and much faster 
than a pre-trained transformer.  

\cite{xlnet19} and \cite{XDHLL19} report 
1.55\% and 1.89\% 
using a pre-trained large transformer, XLnet-large and Bert-large, respectively.  
We observe that the runtime and the amounts of text used for pre-training 
(the last two columns) indicate that their high accuracies come with steep cost at every step: 
pre-training, fine-tuning, and testing.  
Compared with them, an ensemble of GULF-trained DPCNNs is a much lighter-weight solution 
with an appreciable accuracy.  
Also, our ensemble without embedding learning outperforms Bert-base and XLnet-base without pre-training, 
with relatively large differences (Case 2).  
A few attempts of training Bert-large and XLnet-large from scratch also resulted in 
underperforming DPCNNs, 
but we omit the results as we found it infeasible to complete meta-parameter tuning in reasonable time.  

On the other hand, 
it is plausible that the accuracy of the high-performance pre-trained transformers can be further improved 
by applying GULF\ to their fine-tuning, which would further push the state of the art.  
Though currently precluded by our computational constraints, 
this may be worth investigating in the future.  

\subsubsection{Details of the text experiments}
\label{apx:text-details}

\paragraph{Embedding learning}
It was shown in \cite{dpcnn17} that classification accuracy can be improved by 
training an embedding of small text regions (e.g., 3 consecutive words) 
for predicting neighboring text regions (`target regions') on unlabeled data 
(similar to language modeling) 
and then using the learned embedding function to produce additional features for the classifier.  
In this work, 
we trained the following two types of models with respect to use of embedding learning. 
\begin{itemize}
\item Type-0 did not use any additional features from embedding learning. 
\item Type-1 used additional features from the following two types of embedding simultaneously: 
  \begin{itemize}
  \item the embedding of 3-word regions as a function of a bag of words to a 250-dim vector, and 
  \item the embedding of 5-word regions as a function of a bag of word \{1,2,3\}-grams to a 250-dim vector. 
  \end{itemize}
\end{itemize}
Embedding training was done using the entire training set (560K reviews, 391MB) as unlabeled data
disregarding the labels.  

It is worth mentioning that our implementation of embedding learning differs from the 
original DPCNN work \cite{dpcnn17}, as a result of pursuing an efficient implementation in pyTorch 
(the original implementation was in C++). 
The original work used the bag-of-word representation for target regions 
(to be predicted) and minimized squared error with negative sampling. 
In this work we minimized the log loss without sampling where the target probability 
was set by equally distributing the probability mass among the words in the target regions. 

\paragraph{Table \ref{tab:text-results}}
Optimization was done by SGD.  
The learning rate scheduling of the base model and each stage of base-loop\ and 
GULF\ was fixed to 9 epochs with the initial learning rate $\eta$ followed by 
1 epoch with 0.1$\eta$. 
The mini-batch size was 32 for small training data and 128 
for large training data.  
We chose the weight decay parameter from \{1e-4, 2e-4, 5e-4, 1e-3\}
and the initial learning rate from \{0.25, 0.1, 0.05\}, 
using the validation data, 
except that for GULF\ on the large training data, 
we simply used the values chosen for the base model, 
which were weight decay 1e-4 and learning rate 0.1 (with embedding learning)
and 0.25 (without embedding learning).  

For GULF, we chose the number of stages $T$ from \{1,2,\ldots,25\} 
and $\alpha$ from \{0.3, 0.5\}, using the validation data.
$\alpha=0.5$ was chosen in most cases.  

\paragraph{Table \ref{tab:text-results2}}
The ensemble performances were obtained by combining 
\begin{itemize}
\item 
20 DPCNNs ( $T \in$ \{21, 22, \ldots, 25\} $\times$ \{ini:random, ini:base\} $\times$ \{Type-0, Type-1\} )
  in Case 1, and 
\item 
10 DPCNNs ( $T \in$ \{21, 22, \ldots, 25\} $\times$ \{ini:random, ini:base\} $\times$ \{Type-0\}) in Case 2.  
\end{itemize}
To make an ensemble, the model output values were added after softmax.  

\paragraph{Transformers}
The Bert and XLnet experiments were done using HuggingFace's Transformers\footnote{
  {\tt https://huggingface.co/transformers/} 
}
in pyTorch. 
Following the original work, 
optimization was done by Adam with linear decay of learning rate. 
For enabling and speeding up training using a GPU, 
we combined the techniques of gradient accumulation and variable-sized mini-batches (for improving parallelization) 
so that weights were updated after obtaining the gradients from approximately 128 data points. 
128 was chosen, following the original work. 
To measure runtime of transformer testing, we used variable-sized mini-batches for speed-up 
by improving the parallelism on a GPU.  

\section{Introduction}
\label{sec:intro}

This paper presents a new framework to train nonconvex models such as
neural networks. The goal is to learn a vector-valued function 
$\net{\theta}{x}$ that predicts an output $y$ from input $x$,
where $\theta$ is the model parameter. 
For example, 
for $K$-class classification where $y \in \{1,2,\ldots,K\}$, 
$\net{\theta}{x}$ is $K$-dimensional, and it can be linked to 
conditional probabilities via the soft-max logistic function.  
Given a set of training data $S$, %
the standard method for solving this problem is to use 
stochastic gradient descent (SGD) for finding a parameter that minimizes 
on $S$ %
a loss function $L(\net{\theta}{x},y)$ with a regularization term $R(\theta)$: 
$
             \min_\theta \left[ 
                \frac{1}{|S|}\sum_{(x,y)\inS} L(\net{\theta}{x},y)
                + R(\theta)
             \right]. 
$

In this paper, we consider a new framework %
that %
{\em guides training through successive functional gradient descent} so that 
training proceeds with alternating the following: 
\begin{itemize}[itemsep=1pt,topsep=1pt]
\item Generate a guide function so that it is ahead (but not too far ahead) 
      of the current model with respect to the minimization of the loss.  
      This is done by functional gradient descent.   
\item `Push' the model towards the guide function.  
\end{itemize}
Our original motivation was functional gradient learning of additive models in 
{\em gradient boosting} \cite{frie:01}. 
In our framework, essentially, training proceeds with repeating a local search, 
which limits the searched parameter space to the functional neighborhood of the
current parameter at each iteration, instead of searching the entire space at once as the standard method does.   
This is analogous to $\varepsilon$-boosting 
where the use of a very small step-size (for successively expanding the ensemble of weak functions)
is known to achieve better generalization \cite{frie:01}.

For measuring the distances between models, %
we use the Bregman divergence (see e.g., \cite{convex15})
by applying it to the model output.  
Given a convex function $\Bref$, the Bregman divergence $\Bdiv$ is defined by 
\begin{equation}
  \Bdiv(u,v) = \Bref(u) - \Bref(v) - \nabla \Bref(v)^\top (u-v). \label{eq:bdiv}
\end{equation}
This is the difference between $\Bref(u)$ and the approximation of $\Bref(u)$ 
based on the first-order Taylor expansion around $v$.  This means that when $u-v$ is small, 
\begin{equation}
  \Bdiv(u,v) \approx \frac{1}{2} (u-v)^\top(\hessian{\Bref(v)})(u-v) , \label{eq:bdivApprox}
\end{equation}
where $\hessian{h(v)}$ denotes the Hessian matrix of $h$ with respect to $v$. 
Therefore, use of the Bregman divergence has the beneficial effect of utilizing the second-order information. 

We show that the parameter update rule of an induced method 
generalizes that of {\em distillation} \cite{distill14}. 
That is, our framework subsumes iterative {\em self-distillation} 
as a special case.  

Distillation was originally proposed to {\em transfer knowledge} 
from a high-performance but cumbersome model to a more manageable model. 
Various forms of self-distillation, 
which applies distillation to the models of the same architecture, 
has been empirically studied 
\cite{xl19,generation19,LZG18,Ban18,onlinedistill18,deepMutual18,meanTeachers17,gift17}.
One trend is to add to the original scheme, e.g., adding a term to the update rule, 
data distortion/division, more models for mutual learning, and so forth. 
However, we are not aware of any work on theoretical understanding such as
a convergence analysis of the basic self-learning scheme.  

Our theoretical analysis of the proposed framework provides a new functional gradient view of 
self-distillation, and we show a version of the generalized self-distillation
procedure converges to a stationary point of a regularized loss function. 
Our empirical study shows that the
iterative training of the derived method goes through a `smooth path'
in a restricted region
with good generalization performance.  
This is in contrast to standard training, where the entire (and therefore much larger) 
parameter space is directly searched, and thus complexity may not be well controlled.  

\tightpara{Notation}
$\nabla h(v)$ denotes the gradient of a scalar function $h$ with respect to $v$. 
We omit the subscript of $\nabla$ when the gradient is with respect to 
the first argument, e.g., we write $\nabla \net{\theta}{x}$ for $\nabla_\theta \net{\theta}{x}$.  
$\hessian{h(v)}$ denotes the Hessian matrix of a scalar function $h$ with respect to $v$.  
We use $x$ and $y$ for input data and output data, respectively.  
We use $<>$ to indicate the mean, e.g.,
$
  \meanOverXsmall{S}{F(x,y)}
  = \frac{1}{|S|}\sum_{(x,y)\inS} F(x,y). 
$
$L(u,y)$ is a loss function with $y$ being the true output. 
We also let $L_y(u) = L(u,y)$ when convenient. 


\section{Conclusion}

This paper introduces a new framework for guided learning of nonconvex models
through successive functional gradient optimization. A convergence
analysis is established for the proposed approach, and it is shown that
our framework generalizes the popular self-distillation method. Since
the guided learning approach learns nonconvex models in
restricted search spaces, we obtain better generalization performance
than standard training techniques. 

\section*{Acknowledgements}

We thank Professor Cun-Hui Zhang for his support of this research.  


