text:
\section{Conclusion}
A comparison of \ac{arima} and \ac{nnr} for short term solar forecasts showed that it is more likely to find good forecasting models using \ac{nnr}.
Based on this finding, we derived a reduced search space of models that are likely to provide a good \ac{nnr} forecasting model on unseen data.

Future work concerns an extension of the current approach to include exogenous data.
Additionally, other artificial intelligence based methods, e.g., neural networks and support vector regression, are planned to be investigated.

\section{Forecasting Methods}
\label{sec:forecastingMethods}
This section presents the forecasting methods employed in this work.
First, the persistence model used as reference is discussed.
Then, \ac{arima} models and \acl{nnr} are illustrated.

\subsection{Persistence Models}
\label{sec:persistenceModels}
The basic idea of the persistence model is that future values are assumed to be equal to known past values.
In this work, the simple model, $\hat{y}_{t+1|t} = y_t$ is used.
It is also possible to consider a seasonal persistence model using, e.g., data from the previous day, week or year.
This can be described by $\hat{y}_{t+j|t} = y_{t+j-s}$ where $s$ is the seasonal period.

\subsection{Autoregressive Integrated Moving Average (ARIMA) Models}\acused{arima}
\ac{arima} models are widely used for time-series forecasting.
Future values are estimated using a linear combination of previously observed values and forecasting errors.
Also, differencing can be applied to obtain stationary data.
An \ac{arima} model can be described by \cite{BoxJenkinsReinsel1994}
\begin{equation}
	\Phi(B) \nabla^d y_{k} = \theta_0 + \Theta(B) e_{t},
\end{equation}
where
$B$ is the backwards shift operator with $B^m y_{t} = y_{t-m}$,
and
$\nabla$ is the backwards difference with $\nabla^{d} y_{t} = (y_{t}-y_{t-1})^{d}$, \inNLargerEqualZero{d}.
Furthermore,
$\Phi(B) = 1 - \phi_1 B^{1} - \ldots - \phi_p B^{p}$
with $\phi_1, \ldots, \phi_p \in \mathbb{R}$, \inNLargerEqualZero{p}
is the autoregressive part,
and
$\Theta(B) = 1 - \theta_1 B^{1} - \ldots - \theta_q B^{q}$
with $\theta_1, \ldots, \theta_q \in \mathbb{R}$, \inNLargerEqualZero{q}
the moving average part.
Moreover, $e_{t}$ is the difference between measured value and forecast.

To fit a model to data that exhibits seasonality, \ac{arima} models can be extended by seasonality of period $s$.
These, so called \ac{sarima} models, can be described by
\begin{equation}
	\Phi(B) \mathbf{\Phi}(B) \nabla^d \nabla^D_s y_{t} = \theta_0 + \Theta(B) \mathbf{\Theta}(B) e_{t}, \label{eq:sarima}
\end{equation}
where
$\mathbf{\Phi}(B) = 1 - \Phi_1 B^{s} - \ldots- \Phi_P B^{P\,s}$
with $\Phi_1, \ldots, \Phi_P \in \mathbb{R}$, \inNLargerEqualZero{P}
is the seasonal autoregressive part
and
$\mathbf{\Theta}(B) = 1 - \Theta_1 B^{s} - \ldots - \Theta_Q B^{Q\,s}$
with $\Theta_1, \ldots, \Theta_Q \in \mathbb{R}$, \inNLargerEqualZero{Q}
the seasonal moving average part.
Moreover, $\nabla^D_{s} y_{t}= \left(y_{t} -y_{t-s}\right)^D$, \inNLargerEqualZero{D} represents seasonal differencing.

For the implementation of the \ac{arima} models, the econometrics toolbox in MATLAB~2015b was used.
Within the toolbox, maximum likelihood estimation is used to find the model parameters.

\subsection{Nearest Neighbor Regression (NNR)}\acused{nnr}
\label{sec:nnr}

In what follows, based on \cite{Altman1992} we introduce \ac{nnr}.
In \ac{nnr} the pattern of historic data at time instant $t$, $Y_t\in \mathbb{R}^n$, is compared to previous patterns $Y_i \in \mathbb{R}^n$, $i = 0, \ldots, N-1$ in the reference sample
$\mathbb{D} = \left\{ (Y_0, y_1), (Y_1, y_2), \ldots, (Y_{N-1}, y_N) \right\}$,
$N\in\mathbb{N}$
to forecast $\hat{y}_t$.

Although \ac{nnr} allows to use arbitrary elements in $Y_i$, e.g., $Y_i = [y_{i-13} ~ y_{i-111} ~ y_{i-2678}]\tran$, we follow the notion of autoregressive and seasonal autoregressive lags of (seasonal) \ac{arima} models to enable a comparison of both methods.
For \ac{nnr} this means that according to the number autoregressive lags $p$ the vectors $Y_i$ in $\mathbb{D}$ are formed.
Thus, analog to \ac{arima} models the entries of the \ac{nnr} reference sample with $p = 3$ autoregressive lags has the form
$
	Y_{i} = \left[	y_{i} ~ y_{i-1} ~ y_{i-2}\right]\tran.
$

The same holds for the seasonal autoregressive lags.
In a similar fashion as in \eqref{eq:sarima}, the seasonal autoregressive part is multiplied with the autoregressive part.
Consequently, a reference sample with $p = 3$ autoregressive lags, $P = 2$ seasonal autoregressive lags and a season of $s=10$ results in
\[
	Y_{i} = \left[	y_{i} \quad 	y_{i-1} \quad 	y_{i-2} \mid
					y_{i-9} \quad 	y_{i-10} \quad 	y_{i-11} \quad y_{i-12} \mid
					y_{i-19} \quad 	y_{i-20} \quad 	y_{i-21} \quad y_{i-22}\right]\tran.
\]
Thus, as in \eqref{eq:sarima}, the elements of the previous seasons that are associated with $\hat{y}_{i+1|i}$
(here $y_{i-9}$ and $y_{i-19}$) are included in $Y_i$.
This is due to the zero order term, i.e., the $1$, in $\Phi(B)$ and $\mathbf{\Phi}(B)$.

A simple forecasting model can be obtained by combining the \inNLargerZero{k} elements in $\mathbb{D}$ with the smallest distance to $Y_t$.
Using the set of $k$ nearest neighbors of $Y_t$, $\mathbb{K}(Y_t) \subseteq \mathbb{D}$, this can be written as
\begin{equation}
	\hat{y}_{t+1|t} = f(Y_t, \mathbb{D}) = \nicefrac{1}{k}\textstyle\sum_{Y_{i} \in \mathbb{K}(Y_t)} y_{i+1}. \label{eq:nnRegressionFunction}
\end{equation}
Note that the number of neighbors $k$ can be fixed or determined by a maximum distance $\varepsilon \in \mathbb{R}^+$.
With distance $d(Y_i, Y_t)$, the set of neighbors closer than $\varepsilon$ to $Y_t$ is $\mathbb{K}(Y_t) = \{Y_i \in \mathbb{D} \mid d(Y_i, Y_t) \leq \varepsilon \}$.
Note that in this work, the Euclidean norm is used as distance, i.e., $d(Y_i, Y_t) = \|Y_{i} - Y_{t}\|_2$.

The model \eqref{eq:nnRegressionFunction}, can be modified using a weighted average.
In this work, weights inverse to the distance between $x_t$ and the neighbor $x_{i}$, have been considered, i.e.,
\begin{equation}
	\hat{y}_{t+1|t} = f(Y_t, \mathbb{D}) =
	\begin{cases}
		\nicefrac{1}{k}\sum_{Y_{i} \in \mathbb{K}(Y_t)}\frac{1}{d(Y_i, Y_t)} y_{i+1}, & \text{if } d(Y_i, Y_t)\neq 0 ~ \forall Y_i \in \mathbb{K}(Y_t), \\
		y_{l+1}, & \text{if } \exists Y_l \in \mathbb{K}(Y_t) \text{ with } d(Y_t, Y_l)= 0.\\
	\end{cases}
	\label{eq:weightDistance}
\end{equation}

Although, \eqref{eq:nnRegressionFunction} requires no training, the model and the selected reference sample will be referred to as trained model for simplicity.
For implementation, the Statistical Learning Toolbox for MATLAB~\cite{Lin2006} was used.
Note that in this work, approximated nearest neighbors which is faster than the exact nearest neighbors algorithm is used as no significant difference in the forecast accuracy between the two could be observed.

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

In the last years, the capacity of globally installed photovoltaic generators has continuously increased~\cite{Ren21}.
With this growth, the intermittent nature of their infeed has become a major challenge in the operation of electric grids \cite{Denholm2011}.
One important way to address this challenge is to use accurate forecasts to predict the infeed of photovoltaic power plants with a high temporal resolution \cite{West2014}.

Various studies on short term forecasting of
photovoltaic power plant infeed \cite{PedroCoimbra2012, PedroCoimbra2015, Shi2012, Bacher2009, QiaoZeng2012} and
solar irradiance \cite{Reikard2008, SfetsosCoonick2000, Mellit2010} have been published.
Autoregressive integrated moving average (ARIMA)\acused{arima} models have been widely adopted in this context \cite{Reikard2008, PedroCoimbra2012, Bacher2009}.
In the last decades forecasting methods that use techniques from artificial intelligence have become more prominent.
Most of them use \acl{ann} \cite{PedroCoimbra2012, QiaoZeng2012, Reikard2008, SfetsosCoonick2000, Mellit2010}.
Others employ support vector \cite{QiaoZeng2012, Shi2012} and nearest neighbor regression \cite{PedroCoimbra2012, PedroCoimbra2015}.

Despite additional effort that comes with the use of exogenous inputs, e.g., cloud cover or air temperature, only in \cite{PedroCoimbra2012, Bacher2009, SfetsosCoonick2000} forecasts without exogenous inputs are considered.
Moreover, in most publications \cite{PedroCoimbra2012, Bacher2009, Shi2012, SfetsosCoonick2000, Mellit2010, PedroCoimbra2015} the forecasting models are obtained using only one data set.
This makes it hard to draw general conclusions from the analysis that can be transferred to other locations or seasons of the year.
Furthermore, in \cite{PedroCoimbra2012, Mellit2010, Reikard2008, Shi2012, SfetsosCoonick2000, PedroCoimbra2015} the model selection process is not explicitly discussed.
Consequently, even though high forecast accuracies could be achieved for single data sets, due to the missing description of the selection processes, the findings cannot be directly used for different data sets.
For example, in \cite{Mellit2010} the authors state that they used a search and in \cite{Reikard2008} that they tried different model structures to choose an \acl{ann}.
However, they did not provide information on the search space.
In \cite{Shi2012}, a search space for \acl{svr} forecasts is provided but no information on the selection criterion or the final model structure and parameters is given.
The structure of an \acl{ann} forecasting model is examined in \cite{SfetsosCoonick2000} using a sensitivity analysis.
Unfortunately, the publication does not contain sufficient details to replicate the analysis for different data.
To the knowledge of the authors only in \cite{PedroCoimbra2015} a search space is provided and the results are analyzed to gain information about suitable model structures.
However, the analysis is only performed on data from one location.
This makes it hard to draw conclusions that allow identify a reduced search space and in practice often leads to an exhaustive search for suitable models.

In this paper we derive a significantly reduced search space that exhibits a high probability of finding an accurate forecasting model.
As potential forecasting methods, \ac{arima} and \ac{nnr} were considered due to their wide adoption.
Furthermore, their hyperparameters, e.g., the number of autoregressive lags or neighbors, are mostly discrete which allows to form a discrete search space where points in a certain range can be explored.
For each point, i.e., for each set of hyperparameters in this search space, forecasting models are trained and tested using data from different locations and seasons.
Based on the resulting forecast accuracies, conclusions for a reduced search space are drawn.
These include the choice of \ac{nnr} over \ac{arima},
the number of training data points and autoregressive lags
as well as
handling of night data and
use of transmissivity instead of irradiance.
The selection process is explained in much detail allowing others to perform a similar search on different data.

For simplicity and robustness, this work focuses on forecasts of solar irradiance without exogenous inputs using only historic irradiance data, time and location.
Motivated by the smallest intraday interval of energy trading in Germany a \unit[15]{min} prediction step is chosen \cite{EPEX2016}.
To cover a full charge or discharge of a medium size storage unit, a prediction horizon of \unit[3]{h} is~considered.

The remainder of this paper is organized as follows.
Solar irradiance and data preprocessing are discussed in Section~\ref{sec:solarRadiation}.
In Section~\ref{sec:timeSeriesForecasting}, basics on time-series forecasting are discussed.
The forecasting methods considered in this work are presented in Section~\ref{sec:forecastingMethods}.
Then, the model selection procedure is illustrated in Section~\ref{sec:modelSelection}.
Finally, in Section~\ref{sec:analysis} the results of the hyperparameters search are analyzed.

\subsection{Preliminaries}
Throughout this paper, real numbers are denoted by $\mathbb{R}$, nonnegative real numbers by $\mathbb{R}_{0}^+$, positive real numbers by $\mathbb{R}^+$, natural numbers by $\mathbb{N}$ and nonnegative integers by $\mathbb{N}_{0}$.
The transpose of a matrix $x$ is $x^T$ and the Euclidean norm of a vector $x$ is $\| x \|_2$.
The sum over all elements of a set $\mathbb{K}\subset\mathbb{N}$ where each element $i\in\mathbb{K}$ is taken exactly once is denoted $\sum_{i \in \mathbb{K}} x_i$.

\section{Hyperparameter Search}
\label{sec:modelSelection}

\begin{figure}[t]
	\centering

	\includegraphics{figFlowChartGridSearch}%
	\caption{Scheme used to obtain forecast accuracy data for different models.}
	\label{fig:flowChart}
\end{figure}

We aim to compare the forecast accuracy for different forecasting model structures and different sets of historic data.
Therefore, an exhaustive search was performed to evaluate what models approximate the test data best.
In our search, different model structures are trained using different sets of historic data.
To identify model structure and historic data, hyperparameters are used.
Moving along Fig.~\ref{fig:flowChart}, in this section first training and test data are discussed.
Then, the hyperparameter search space is posed.
Finally, training and test of the models is illustrated.

\subsection{Training and Test Data}
\label{sec:data}

For training and test data sets \cite{NRELNevada,NRELSRRL,ARM} were used.
From each data set three different seasons were selected to include different climatic situations in the hyperparameter search.
For every location and every season, i.e., for 9 data sets, the observations were divided into training data (up to 2~months) and independent test data (1~week) that was not used for training.
The test data was chosen to always follow directly the data used for training.
As shown in Fig.~\ref{fig:testdata}, the nine weeks of test data include different climatic situations, e.g.,
a sunny week in December from the \ac{nrel} Clark Station,
a cloudy week in March from the \ac{arm} Facility, and
various weeks with sunny, foggy and cloudy days.

Note that all data sets \cite{NRELNevada,NRELSRRL,ARM} have a temporal resolution of \unit[1]{min}.
The \unit[15]{min} average values $y_1, \ldots, y_N$ were derived from the original time series $y'_1, \ldots, y'_{15N}$ as $y_t = \nicefrac{1}{15} \textstyle\sum_{q=1}^{15} y'_{q+15(t-1)}$.

\begin{figure}[h]

	\includegraphics{dataTest}%
	\caption[Test Data]{Measured solar irradiance data used to test the trained models.}
	\label{fig:testdata}
\end{figure}

\subsection{Hyperparameters}
\label{sec:hyperparameters}

\subsubsection{Data Hyperparameters.}
\label{sec:dataHyperparameters}

The hyperparameters that concern the training data are independent of the forecasting method and the model structure.
We considered modifications of training data by
\begin{enumerate*}
	\item pre- and postprocessing,
	\item handling of night data ans
	\item number of data points used.
\end{enumerate*}

1. Regarding the pre- and postprocessing, two approaches were considered.
(a) The time series of irradiance, $I_t$, is directly used to train the models without further preprocessing.
(b) Transmissivity $\tau_t = \nicefrac{I_t}{I_{\text{e},t}}$ is derived from $I_t$ and extraterrestrial irradiance~$I_{\text{e},t}$.
Then, the models are trained to forecast $\hat{\tau}_{t+j|t}$.
Using $\hat{\tau}_{t+j|t}$ the irradiance forecast is then calculated via $\hat{I}_{t+j|t} = I_{\text{e},t} \, \hat{\tau}_{t+j|t}$.

2. Three cases related to the handling of data points at night were considered:
(a) include all data points from day and night,
(b) only include data points between \unit[5]{am} and \unit[8]{pm} local time, and
(c) only include data points between sunrise and sunset, i.e., observations with $\zeta_t > 90.83^\circ$ (see Section~\ref{sec:handlingOfNightData}).
As it can be cumbersome to implement models with seasonality for a varying number of data points per day, case (3) was not considered for seasonal models.

3. Regarding the number of training data points, $1$, $3$, $7$, $14$, $30$, and $60$ days of data were considered.
For seasonal models, the data was chosen to include at least $14$ days.

\subsubsection{Model Hyperparameters.}

\input{tableParameterSpaceCompact.tex}

For every forecasting method, the model structure can be uniquely described using model specific hyperparameters.
Note the difference to model parameters of a trained model.
For example, the number of autoregressive lags $p$ is a hyperparameter.
The coefficients of the autoregressive part, $\phi_1, \ldots, \phi_p$, of a trained \ac{arima} model are model parameters.
The subsets of hyperparameters considered in the search for are shown in Table~\ref{tab:parameter}.
Note that in Tables~\ref{tab:parameter}(\subref{tab:parameterNnr}) and~\ref{tab:parameter}(\subref{tab:parameterSeasonalNnr}) either the number of neighbors $k$ or the maximal distance to the neighbors $\varepsilon$ is used.
Further note that the season $s$ was chosen to be \unit[1]{d} for \ac{sarima} and \ac{snnr}.

\subsection{Model Training}
For each combination of hyperparameters described in Section~\ref{sec:hyperparameters}, we trained nine models, i.e., one for each location and each season (see Section~\ref{sec:data}).
This lead to more than 250,000 trained forecasting models with different hyperparameters.
Note that some combinations of hyperparameters did not result in suitable models.
Especially for \ac{arima}, no stable model could be derived for certain hyperparameters.
Furthermore, for \ac{sarima} \unit[3.3]{\%} of the hyperparameters were excluded, mostly due to very long training procedures that rendered them unsuitable for practical use.

\subsection{Model Test}
\label{sec:modelTest}

In the model test, the performance of the trained models was evaluated.
Therefore, every trained model was tested by forecasting irradiance of the unknown test data (see Fig.~\ref{fig:testdata}) that followed the training data.
Due to zero irradiance, data points during the night were excluded from the evaluation based on the zenith angle $\zeta_t$ (see Section~\ref{sec:handlingOfNightData}).
For each of the data point in the test data, a $12$ step ahead, i.e., \unit[3]{h}, prediction was obtained.
Using the forecast values and the test data, the $\operatorname{RMSE}_j$ (see Section~\ref{sec:rmse}) was then calculated for $j=1, \ldots, 12$ and stored.
In the next section, we will analyze the resulting $\operatorname{RMSE}_j$ for every combination of hyperparameters, location and season to identify models that are likely to provide good forecasts.

\section{Analysis}
\label{sec:analysis}
In this section, we analyze the results of the hyperparameter search in Section~\ref{sec:modelSelection}.
Our goal is to draw conclusions that help to reduce the space of future hyperparameter searches.
More precisely, we aim to find forecasting methods and ranges of hyperparameters that increase the probability of finding a suitable model for unseen data.
Before starting the analysis, some remarks are posed.

\begin{remark}\label{rem:referenceModel}
	As reference, we use the persistence model on transmissivity as it outperformed the persistence model on irradiance for all test data sets.
	In the reference model historic transmissivity is derived using $\tau_t = \nicefrac{I_t}{I_{\text{e},t}}$.
	Then a persistence forecast is performed, i.e., we set $\hat{\tau}_{t+j|t} = \tau_t$ for $j = 1, \ldots, J$.
	Finally, the irradiance forecast is determined as $\hat{I}_{t+j|t} = \hat{\tau}_{t+j|t} \, I_{\text{e},t+j}$ for $j = 1, \ldots, J$.
\end{remark}

\begin{remark}
	The style of the box plots used throughout this work is illustrated in in Fig.~\ref{fig:boxPlotExplaination}.
	Here, $m$ marks the median.
	The box around the median contains all data from the 25th ($q_1$) to the 75th ($q_3$) percentile.
	The left whisker marks the lowest occurring value within $q_1 - 1.5 (q_3 - q_1)$ and the right whisker marks the highest occurring value within $q_3 + 1.5 (q_3 - q_1)$.
	Due to numerous outliers, only every 100th outlier is shown.
	To assure that the forecast with the highest accuracy is considered the lowest outlier is always included.

	\begin{figure}[h]
		\centering

		\includegraphics{boxPlotExplaination}%
		\caption[/boxPlotExplaination.tex]{Box plot with median $m$, quartiles $q_1$ and $q_3$, whiskers and outliers.}
		\label{fig:boxPlotExplaination}
	\end{figure}
\end{remark}

\begin{remark}
	Most distributions of the forecast accuracy discussed in this section are similar for the different test data sets shown in Fig.~\ref{fig:testdata}.
	Therefore, most distributions were combined in one plot that includes the values of all locations.
\end{remark}

\begin{remark}
	Throughout this paper, the forecast at prediction steps $j = 1, 4, 12$ are analyzed to approximate the evolution of the forecast accuracy over prediction horizon $J=12$.
	For the \unit[15]{min} sampling time considered, they correspond to a forecast \unit[15]{min}, \unit[1]{h} and \unit[3]{h} ahead.
\end{remark}

\subsection{Forecasting Method}
\label{sec:forecastMethod}

We first compare the forecasting methods \ac{arima} and \ac{nnr} along with the persistence model on transmissivity that acts as a reference (see Remark~\ref{rem:referenceModel}).
In Fig.~\ref{fig:boxPlotArimaNnr}, box plots illustrating the distributions of forecast accuracies of all models included in the search are shown.

\begin{figure}[htbp]
	\centering

	\includegraphics{boxPlotArimaNnr}%
	\caption[/boxPlotArimaNnr.tex]{\ac{rmse} for \unit[15]{min}, \unit[1]{h} and \unit[3]{h} ahead predictions with \ac{arima} and \ac{nnr} models.}
	\label{fig:boxPlotArimaNnr}
\end{figure}

The box plots in Fig.~\ref{fig:boxPlotArimaNnr} show that for predictions \unit[15]{min} ahead, almost no model outperforms the persistence model.
For a prediction step of \unit[1]{h} some models outperform the persistence forecaster.
Furthermore, the \ac{rmse} of forecasts performed with \ac{nnr} varies less than the \ac{rmse} of \ac{arima} forecasts.
This indicates a higher probability of finding a good forecasting model using \ac{nnr}.
For predictions of \unit[3]{h} ahead, the potential improvement over the persistence model is significant for all methods.
As already observed for \unit[1]{h} ahead predictions, the \ac{rmse} of \ac{nnr} varies less.
Additionally, the 25th percentile is remarkably lower for \ac{nnr}.
The highest forecast accuracy achieved is very alike for \ac{nnr} and \ac{arima} for all prediction steps.
Still, \ac{nnr} achieves much lower medians and lower 25th and 75th percentiles.
Consequently, the probability of finding a suitable \ac{nnr} model is much higher.
Therefore, \ac{arima} and \ac{sarima} models will be disregarded in the next steps of the analysis.
Consequently, in what follows we only investigate how the search spaces for \ac{nnr} and \ac{snnr} can be reduced to increase the probability of finding an accurate forecasting model.

\subsection[Preprocessing]{Preprocessing of \ac{nnr} and \ac{snnr}}
\label{sec:preprocessing}

As stated in Section~\ref{sec:dataHyperparameters}, models using irradiance or the normalized transmissivity values were considered in the search.
In Fig.~\ref{fig:boxPlotGhiCsi} the forecast accuracy of both approaches is shown.
It can be observed that the box plots for \unit[15]{min} forecasts are very alike.
Also, little difference can be seen for \unit[1]{h} ahead predictions.
However, for \unit[3]{h} ahead predictions the transmissivity based forecasts often outperform the irradiance based forecasts.
The only exception from this observation could be found for data from the \ac{nrel} Clark Station during March 2--8 (see Fig.~\ref{fig:boxPlotNnrPreprocessingException}), where the best models use irradiance data.
Still the lower whiskers and the medians in this case are very alike.

\begin{figure}[ht]
	\centering

	\includegraphics{boxPlotNnrPreprocessing}%
	\caption[Box plot nnr preprocessing]{%
	\ac{rmse} of \ac{nnr} and \ac{snnr} using different preprocessing for all locations and seasons.}
	\label{fig:boxPlotGhiCsi}
\end{figure}

The box plots indicate that using transmissivity yields a higher probability of finding an accurate model.
Although for data set \cite{NRELNevada} a better irradiance model for the \unit[3]{h} prediction could be found, the overall distribution in this case shows no significant additional difference to Fig.~\ref{fig:boxPlotGhiCsi}.
As the single best model is hard to find, we are confident that using transmissivity yields a higher probability of finding a good model.
Therefore, irradiance forecasts are excluded from the reduced search space.

\begin{figure}[h!]
	\centering

	\includegraphics{boxPlotNnrPreprocessingException}%
	\caption[Box plot nnr preprocessing exception]{\ac{rmse} of \ac{nnr} and \ac{snnr} using different preprocessing at \ac{nrel} Clark Station, March 2--8.}
	\label{fig:boxPlotNnrPreprocessingException}
\end{figure}

\subsection[Handling of Night Data]{Handling of Night Data of \ac{nnr} and \ac{snnr}}
\label{sec:nightData}

As stated in Section~\ref{sec:modelTest}, data points during the night were not included in the evaluation of the forecast accuracy.
However, in some cases they were used in the data of the reference sample to allow for smooth transitions between two days.
As stated in Section~\ref{sec:data}, we considered training data that
(a) includes the data of the entire day (and night),
(b) only includes data points from sunrise to sunset and
(c) only includes data points from \unit[5]{am} to \unit[8]{pm}.

\begin{figure}[t]
	\centering

	\includegraphics{boxPlotNnrHandlingNightData}%
    \caption[Box Plot NNR Night Data]{%
	\ac{rmse} of \ac{nnr} and \ac{snnr} for different handling of night data: (a)~entire day, (b)~sunrise to sunset and (c)~\unit[5]{am} to \unit[8]{pm}.}
	\label{fig:boxPlotNightData}
\end{figure}

In Fig.~\ref{fig:boxPlotNightData}, the forecast accuracy for different handling of night data is shown.
It can be observed that the forecast accuracy is similar for the \unit[15]{min} forecast.
For forecasts \unit[1]{h} ahead, the difference between the approaches remains small.
However, for forecasts \unit[3]{h} ahead, using all data points for training, i.e., including night values, outperforms the other approaches.
Therefore, all data points including night data are considered in the reduced search space.

\subsection[Size of the Reference Sample]{Size of the Reference Sample of \ac{nnr} and \ac{snnr}}
\label{sec:sizeSample}

The reference sample is chosen such that it includes that last values of historic data up to the most recent data point.
A large reference sample, i.e., a set $\mathbb{D}$ with a high number of elements, therefore includes data that reaches further into that past than a small reference sample.

\begin{figure}[ht]
	\centering

	\includegraphics{fanChartNnrSizeDataTraining}%
	\caption[Fan chart NNR size data training]{%
	\ac{rmse} of \ac{nnr} and \ac{snnr} using different sizes of the reference sample.}
	\label{fig:sizeDataTraining}
\end{figure}

In Fig.~\ref{fig:sizeDataTraining}, the forecast accuracy using a reference sample of \unit[1]{d}, \unit[3]{d}, and \unit[7]{d} as well as two weeks (\unit[14]{d}), three weeks (\unit[21]{d}), one month (\unit[30]{d}) and two months (\unit[60]{d}) is shown.
The plots for all prediction steps shows a small decrease of the 25th and the 75th percentile as well as the lowest \ac{rmse} for an increasing number of elements in the reference sample.
Considering the largest \ac{rmse} of each prediction step, the inverse effect can be observed.
In general little difference can be observed in the forecast accuracy, when varying the size of the reference sample.
However, as the forecast accuracy increases using more data, a reference sample of \unit[60]{d} is used in the reduced search space.

\subsection[Autoregressive Lags]{Autoregressive (AR) Lags of \ac{nnr} and \ac{snnr}} \acused{ar}
\label{sec:analysisLags}
This section focuses on the \ac{ar} lags, i.e., the data points of the historic data used in the elements of the reference sample.
We differentiate between \ac{ar} and seasonal \ac{ar} lags.
The former refer to \unit[15]{min} steps prior to the first forecast time instant.
The latter refer to \unit[24]{h} steps, used to include the days prior to the forecast.
\begin{figure}[t]
	\centering

	\includegraphics{fanChartNnrLags}%
	\caption[fanChartNnrLags]{\ac{rmse} of \ac{nnr} using different \ac{ar} lags.}
	\label{fig:fanChartNnrLags}
\end{figure}
It can be observed in Fig.~\ref{fig:fanChartNnrLags} that for all prediction steps the models with the highest forecast accuracy has a small number of \ac{ar} lags.
However, no clear tendency can be identified that supports choosing a particular number of \ac{ar} lags.

\begin{figure}[ht]
	\centering

	\includegraphics{fanChartNnrSeasonalLags}%
	\caption[fanChartNnrSeasonalLags]{%
	\ac{rmse} of \ac{snnr} using different seasonal \ac{ar} lags.
	For every seasonal \ac{ar} lag, the non-seasonal \ac{ar} lags increase from $1$ on the left side of the plot to $11$ on the right side of the small subplot.}
	\label{fig:fanChartNnrSeasonalLags}
\end{figure}

Fig.~\ref{fig:fanChartNnrSeasonalLags} shows the \ac{rmse} for different numbers of seasonal \ac{ar} lags.
Here, for every seasonal \ac{ar} lag, the number of \ac{ar} lags was increased from $1$ (left) to $11$ (right).
It can be observed that for \unit[15]{min} forecasts the \ac{rmse} slightly increases with the number of seasonal \ac{ar} lags.
Furthermore, for every seasonal \ac{ar} lag, the \ac{rmse} also increases with increasing number of non-seasonal \ac{ar} lags.
In contrast, for the \unit[3]{h} forecasts, the \ac{rmse} slightly decreases with the number of seasonal lags.
Comparing Figs.~\ref{fig:fanChartNnrLags} and~\ref{fig:fanChartNnrSeasonalLags}, shows that the variance decreases for the \unit[1]{h} and \unit[3]{h} forecasts with increasing number of seasonal lags.
This results in a higher probability of finding a good model for a higher number of seasonal \ac{ar} lags.
Unfortunately, the \unit[15]{min} forecasts do not show the same effect.
Still, as the \unit[15]{min} forecast barely outperform the persistence model it seems reasonable focus on the improvement for larger horizons and consider only seasonal models in the reduced search space.

\subsection[NNR Weights]{Weights of \ac{nnr} and \ac{snnr}}
\label{sec:weights}
As stated in Section~\ref{sec:hyperparameters} uniform weights and weights inverse to the distance to the neighbors were considered.
In the following, we analyze the effect of the distance to neighbors on the distribution of the \ac{rmse}.
\begin{figure}[bt]
	\centering

	\includegraphics{boxPlotNnrUniformDistance}%
	\caption[Boxplot weights knn]{%
	\ac{rmse} using uniform weights and weights inverse to the distance of the neighbor.}
	\label{fig:boxPlotKnnWeights}
\end{figure}
As can be seen in Fig.~\ref{fig:boxPlotKnnWeights}, the box plots hardly differ for all prediction steps.
This indicates that there is no strong correlation between forecast accuracy and weights.
However, because of an easier implementation, uniform weights were chosen for the reduced search space.

\subsection[Definition of the Neighborhood]{Definition of the Neighborhood of \ac{nnr} and \ac{snnr}}
\label{sec:neighborhood}
As stated in Section~\ref{sec:nnr}, the neighborhood can be chosen to be a fixed number of neighbors, $k$, or a maximum distance, $\varepsilon$.
\begin{figure}[t]
	\centering

	\includegraphics{boxPlotKnnKernel}%
	\caption[BoxplotKnnKernel]{%
	\ac{rmse} using a threshold $\varepsilon$ or a fixed number of neighbors $k$ to define the neighborhood $\mathbb{K}$.}
	\label{fig:boxPlotKnnVsKernel}
\end{figure}
The comparison in Fig.~\ref{fig:boxPlotKnnVsKernel} shows that the highest accuracy achieved with both hyperparameters is very alike for all prediction steps.
However, the 25th and 75th percentile, the upper whisker and the median are lower for models using a fixed number of neighbors $k$.
Therefore, in the reduced search space, only models using a fixed number of neighbors $k$ are considered.

\subsection[Number of Neighbors]{Number of Neighbors $k$ of \ac{nnr} and \ac{snnr}}
\label{sec:numberNeighbors}

In Fig.~\ref{fig:fanChartNumberOfNeighbors}, the forecast accuracy of all \acl{knn} models included in the search are shown.
It can be observed, that the forecast accuracy mostly changes using $1$ to $5$ neighbors.
The distribution of the \ac{rmse} becomes narrower for most cases with an increasing $k$.
Furthermore, the lowest \ac{rmse} decreases slightly for $k \geq 10$.
Therefore, models using between 10 and 20 neighbors are considered the reduced search space.

\begin{figure}[htb]
	\centering

	\includegraphics{fanChartNnrNumberOfNeighbors}%
	\caption[Fan chart kNN - number of neighbours k.]{\ac{rmse} of \ac{nnr} and \ac{snnr} using a different number of neighbors $k$.}
	\label{fig:fanChartNumberOfNeighbors}
\end{figure}

\subsection{Summary}
\label{sec:results:summary}

In summary, the following conclusions can be drawn from the analysis.
\begin{enumerate}
	\item The probability of finding a sufficiently accurate model to forecast irradiance is much higher for \ac{nnr} and \ac{snnr} than for \ac{arima} and \ac{sarima} models (see Section~\ref{sec:forecastMethod}).
	\item Data pre- and postprocessing of \ac{nnr} and \ac{snnr} models:
	\begin{enumerate}
		\item Use transmissivity instead of irradiance (Section~\ref{sec:preprocessing}).
		\item Include night data in the reference sample (Section~\ref{sec:nightData}).
		\item Use data from the past \unit[60]{d} in the reference sample (Section~\ref{sec:sizeSample}).
	\end{enumerate}
	\begin{samepage}
	\item Hyperparameters of \ac{nnr} and \ac{snnr} models: \label{enum:results:hyperparameters}
	\begin{enumerate}
		\item Seasonal models are beneficial (Section~\ref{sec:analysisLags}).
		\item Favor uniform weights over weights inverse to distance (Section~\ref{sec:weights}).
		\item Choose the neighbors using a fixed number $k$ instead of a maximum distance (Section~\ref{sec:neighborhood}).
		\item Search for models with 10 to 20 neighbors (Section~\ref{sec:numberNeighbors}).
	\end{enumerate}
	\end{samepage}
\end{enumerate}

Based on these findings the hyperparameter search space, i.e., the number of potential models could be reduced from more than 250,000 to less than 1,000.
The search space for suitable \ac{knn} forecasters now only includes the
number of autoregressive lags with range $1, 2, \ldots, 11$,
number of seasonal autoregressive lags with range $1, 2, \ldots, 7$ and
number of neighbors with range~$10, 11, \ldots, 20$.

\section{Solar Irradiance}
\label{sec:solarRadiation}
The sunlight reaching the outer earth's atmosphere is called extraterrestrial solar radiation.
It can be estimated from the energy emitted by and the position of the sun.
The sunlight reaching a horizontal plane on earth per unit area at time $t$, is called global horizontal irradiance $I_t\in\R_{0}^+$.
It includes direct normal irradiance, which originates directly from the sun, and diffuse irradiance, which includes scattered and ground reflected components, \cite{ISO9488}.

\begin{figure}[tbp]
	\centering

	\includegraphics{figSolarIrradiance}%
	\caption[Irradiance, Transmissivity]{Global horizontal irradiance, extraterrestrial irradiance and transmissivity over a duration of \unit[5]{d}.
	The irradiance was taken from \cite{ARM}, the extraterrestrial irradiance was estimated using \cite{RA2004}.
	 \label{fig:solarIrradiance}}
\end{figure}

For forecasting, it can be beneficial to normalize the solar irradiance $I_t$ by the extraterrestrial irradiance $I_{\text{e},t}\in\R_{0}^+$ (see Fig.~\ref{fig:solarIrradiance}).
This yields to transmissivity defined as
$
	\tau_t = \nicefrac{{I_t}}{{I}_{\text{e},t}}
$
\cite{QiaoZeng2012}.
In this work, the solar position algorithm from \cite{RA2004} is used to estimate the extraterrestrial irradiance $I_\text{e}(t)$.
It only requires the position of the respective surface to calculate the zenith angle $\zeta_t\in[0, 2\pi)$, i.e., the incidence angle of the sun light on a horizontal plane on the earth's surface \cite{Quaschning2009}.
Using $\zeta_t$, the extraterrestrial irradiance can be determined by
$
	I_{\text{e},t} = \epsilon_t I_{\text{s}} \cos(\zeta_t),
$
where $\epsilon_t$ is a correction factor and $I_{\text{s}} = \unitfrac[1360.8]{W}{m^2}$ is the solar constant.

\section{Time-Series Forecasting}
\label{sec:timeSeriesForecasting}

A time series is a collection of \inNLargerZero{n} chronologically ordered observations $y_1, y_2, \ldots, y_n$.
In this collection, every element $y_t$ refers to an observation performed at a time instant $t = 1, \ldots, n$.
In this work, we focus on univariate forecasts, i.e., forecasts using only historic irradiance to forecast future irradiance.
Thus, all elements $y_t\in\mathbb{R}$ for $t = 1, \ldots, n$ are real-valued scalars.

Broadly speaking, a forecasting method is a procedure to estimate future values $\hat{y}_{t+j|t}$ of a time series based on past values $Y_t \in \mathbb{R}^n$, i.e.,
$
	\hat{y}_{t+j|t} = f(Y_t).
$
Here, $\hat{y}_{t+j|t}$ refers to the value at prediction step $\inNLargerZero{j}$ with last known value at time $t\in\mathbb{N}$.
In case of univariate forecasting, $Y_t$ only includes the present and past values of the time series being forecast, i.e.,
$
	Y_t = [y_{t} ~ y_{t-1} ~ \cdots ~ y_{t-n+1}]\tran
$.

The vector $Y_t$ can also be based on non-consecutive lags, e.g., to represent seasonal behavior.
An example for such a vector is
$
	Y_t = [y_{t} ~ y_{t-s+1} ~ y_{t-s}]\tran,
$
where \inNLargerZero{s} is the length of the season, e.g., one day.
In this example,
the present observation, $y_t$,
the observation from the last day, $y_{t-s+1}$, at the same time as the predicted value, $\hat{y}_{t+1}$,
and from the current time minus \unit[24]{h}, $y_{t-s}$, are included.

\subsection{Multi-Step Forecasts}
\label{sec:multipleStepForecasts}
\begin{subequations}
To create multi-step ahead forecasts, a recursive strategy can be applied \cite[Section 4.]{Bontempi2013}.
The forecasting model is repeatedly used for one step forecasts, where the elements of $Y_t$ are adapted at each step until prediction step \inNLargerZero{J} is reached.
For example, if $Y_t$ includes the last three data points, $Y_t = [y_{t} ~ y_{t-1} ~ y_{t-2}]^{T}$
then the procedure for $J = 12$ prediction steps would be
\begin{align}
	\hat{y}_{t+1|t}		& = f(Y_t) = f([y_{t} \quad y_{t-1} \quad y_{t-2}]\tran), \\
	\hat{y}_{t+2|t}		& = f(\widehat{Y}_{t+1|t}) = f([\hat{y}_{t+1|t} \quad y_{t} \quad y_{t-1}]\tran), \\
						& \mathrel{\makebox[\widthof{=}]{\vdots}} \nonumber\\
	\hat{y}_{t+12|t}	& = f(\widehat{Y}_{t+11|t}) = f([\hat{y}_{t+11|t} \quad \hat{y}_{t+10|t} \quad \hat{y}_{t+9|t}]\tran).
\end{align}
For seasonal models, the vector  $\widehat{Y}_{t+j|t} \in \mathbb{R}^n$ is equally adapted at each step in a similar manner.
\end{subequations}

\subsection{Evaluation of Forecast Accuracy}
\label{sec:evaluation}

\subsubsection{Handling of Night Data.}
\label{sec:handlingOfNightData}

Forecasts of data points during the night, i.e., of zero irradiance values, were not included in the evaluation of the forecast accuracy.
Therefore, observations with a zenith angle
$
	\zeta_t \leq 90.83^\circ
$
were excluded based on \cite{Meeus1991}.
The angle $\zeta_t$ was taken from the data sets if available.
Otherwise it was estimated using \cite{RA2004}.

\subsubsection{Root Mean Square Error (RMSE).}\acused{rmse}
\label{sec:rmse}
For each prediction step $j = 1,\ldots, J$ over prediction horizon $J$, the \ac{rmse} is calculated to analyze the error of each step separately.
Assuming that \inNLargerZero{m} forecasts were performed, the \ac{rmse} of prediction step $j$ is
\begin{equation}
	\label{eq:rmse}
	\operatorname{RMSE}_j = \textstyle\sqrt{\frac{1}{m} \sum_{t=1}^{m} \left(\hat{y}_{t+j|t}-y_{t+j}\right)^2}.
\end{equation}


text:
\section{Supplementary Materials for Online hyper-parameter optimization by real-time recurrent learning}

\iffalse
\begin{figure*}[H]
    \centering
    \begin{minipage}{0.49\textwidth}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLr_wrt_lr_curve_comp_sgd_mlr0.000050.pdf}
    \subcaption*{Mlr=5e-5}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLr_wrt_lr_curve_comp_sgd_mlr0.000010.pdf}
    \subcaption*{Mlr=1e-5}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLr_wrt_lr_curve_comp_sgd_mlr0.000005.pdf}
    \subcaption*{Mlr=5e-6}
    \end{minipage}
    \subcaption{Influence matrix norm w.r.t $\alpha$}
    \end{minipage}
    \begin{minipage}{0.49\textwidth}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLl2_wrt_l2_curve_comp_sgd_mlr0.000050.pdf}
    \subcaption*{Mlr=5e-5}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLl2_wrt_l2_curve_comp_sgd_mlr0.000010.pdf}
    \subcaption*{Mlr=1e-5}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_dPdLl2_wrt_l2_curve_comp_sgd_mlr0.000005.pdf}
    \subcaption*{Mlr=5e-6}
    \end{minipage}
    \subcaption{Influence matrix norm w.r.t $\lambda$}
    \end{minipage}
    \caption{L2 norm of influence matrix w.r.t learning rate $\alpha$ and weight decay coefficient $\lambda$ VS meta Learning rate}
    \label{fig:normJ_fix_mlr}

    \begin{minipage}{\textwidth}
    \begin{minipage}{0.245\textwidth}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_dPdLr_wrt_mlr_curve_comp_sgd_lr0.001000.pdf}
        \subcaption*{ILr=0.0001}
        \end{minipage}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_dPdLr_wrt_mlr_curve_comp_sgd_lr0.000100.pdf}
        \subcaption*{ILr=0.001}
        \end{minipage}
    \subcaption{Influence matrix norm w.r.t $\alpha$}
    \end{minipage}
    \begin{minipage}{0.245\textwidth}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_dPdLl2_wrt_mlr_curve_comp_sgd_lr0.001000.pdf}
        \subcaption*{ILr=0.0001}
        \end{minipage}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_dPdLl2_wrt_mlr_curve_comp_sgd_lr0.000100.pdf}
        \subcaption*{ILr=0.001}
        \end{minipage}
    \subcaption{Influence matrix norm w.r.t $\lambda$}
    \end{minipage}
    \begin{minipage}{0.245\textwidth}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_mlr_dPdLr_curve_comp_sgd_lr0.001000.pdf}
        \subcaption*{ILr=0.0001}
        \end{minipage}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_mlr_dPdLr_curve_comp_sgd_lr0.000100.pdf}
        \subcaption*{ILr=0.001}
        \end{minipage}
        \subcaption{Influence norm w.r.t $\alpha$}
    \end{minipage}
    \begin{minipage}{0.245\textwidth}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_mlr_dPdLl2_curve_comp_sgd_lr0.001000.pdf}
        \subcaption*{ILr=0.0001}
        \end{minipage}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_mlr_dPdLl2_curve_comp_sgd_lr0.000100.pdf}
        \subcaption*{ILr=0.001}
        \end{minipage}
    \subcaption{Influence norm w.r.t $\lambda$}
    \end{minipage}
    \end{minipage}
    \caption{L2 norm of Influence w.r.t learning rate $\alpha$ and weight decay coefficient $\lambda$ VS meta Learning rate}
    \label{fig:normJ_lr}
\end{figure*}
\fi

\section{Background}

\subsection{Iterative optimization}

A standard technique for training large-scale machine learning models, such as deep neural networks (DNN), is stochastic optimization. In general, this takes the form of parameters changing as a function of the the data and the previous parameters, iteratively updating until convergence:
\begin{align}

    \boldsymbol{\theta}_{\tau+1} = g(\boldsymbol{\theta}_{\tau}, \mathbf{B}_{\tau}; \boldsymbol{\varphi}),
    \label{eq:update_rule}
\end{align}
where $\boldsymbol{\theta}$ denotes the parameters of the model, $g$ is the update rule reflecting the learning objectives, ${\boldsymbol \varphi}$ denotes the hyperparameters,  and $\mathbf{B}_{\tau}$ is the random subset of training data (the minibatch) at iteration $\tau$.

As a representative example, when training a DNN $f$ using stochastic gradient descent, (SGD) learning follows the loss gradient approximated using a random subset of the training data at each time step
\begin{align}
\label{eq:sgd}
    \boldsymbol{\theta}_{\tau+1} \leftarrow \boldsymbol{\theta}_{\tau} - \frac{\alpha}{|\mathbf{B}_\tau|} \sum_{(x,y) \in \mathbf{B}_\tau} \nabla_{\theta} l(y, f(x; \boldsymbol{\theta}_{\tau})),
\end{align}
where $l(\cdot)$ is the per-example loss function and $\alpha$ is the learning rate.
In this case, $g(\cdot)$ is a simple linear update based on the last parameter and the scaled gradient. Other well-known optimizers, such as  RMSprop~\citep{Tieleman2012} or Adam~\citep{Kingma2015}, can also be expressed in this general form.

\subsection{Recurrent neural networks}
\label{sec:rtrl}

An RNN is a time series model. At each time step $t =\{ 1\cdots T\}$, it updates the memory state $\mathbf{h}_t$ and generates the output $\mathbf{o}_{t+1}$ given the input $\mathbf{x}_t$:
\begin{align}
\label{eq:rnn}
   \left( \mathbf{h}_{t+1}, \mathbf{o}_{t+1} \right) = r(\mathbf{h}_{t}, \mathbf{x}_t; \boldsymbol{\phi}),
\end{align}
where the recurrent function $r$ is parametrized by $\boldsymbol{\phi}$.

Learning this model involves optimizing the total loss over $T$ steps with respect to $\boldsymbol{\phi}$ by gradient descent. BPTT is typically used to calculate the gradient by unrolling the network dynamics and backward differentiating through the chain rule:
\begin{align*}
    \nabla_\phi \mathcal{L}

                            &= \sum^T_{t=1}  \Big(\sum^T_{s\geq t+1}\frac{\partial \mathcal{L}_{s}}{\partial \mathbf{h}_{t+1}} \frac{\partial \mathbf{h}_{t+1}}{\partial \mathbf{h}_{t}} + \frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_{t}} \Big)
                            \frac{\partial \mathbf{h}_t}{\partial \boldsymbol{\phi}},
\end{align*}
where $\mathcal{L}_t$ is the instantaneous loss at time $t$, and the double summation indexes all losses and all the applications of the parameters.
Because the time complexity for BPTT grows with respect to $T$, it can be challenging to compute the gradient using BPTT when the temporal horizon is lengthy. In practice, this is mitigated by truncating the computational graph (`truncated BPTT').

\subsection{Real-time recurrent learning}

RTRL \citep{Williams1989} is an online alternative to BPTT, which --instead of rolling out the network dynamics-- stores a set of summary statistics of the network dynamics that are themselves updated online as new inputs come in. To keep updates causal, it uses a forward view of the derivatives, instead of backward differentiation, with the gradient computed as
\begin{align}
    \nabla_{\boldsymbol{\phi}} \mathcal{L} &

     = \sum^T_{t=1}  \frac{\partial \mathcal{L}_{t+1}}{\partial \mathbf{h}_{t+1}} \Big(\frac{\partial \mathbf{h}_{t+1}}{\partial \mathbf{h}_{t}}\frac{\partial \mathbf{h}_{t}}{\partial \boldsymbol{\phi}} +\frac{\partial \mathbf{h}_{t+1}}{\partial \boldsymbol{\phi}_t} \Big).
     \label{eq:rtrl}
\end{align}
The Jacobian $\Gamma_{\tau}=\frac{\partial \mathbf{h}_{t}}{\partial \boldsymbol{\phi}}$ (also referred to as the {\em influence matrix}), dynamically updates as
\begin{align*}
    \Gamma_{\tau} = D_\tau \Gamma_{\tau-1} + G_{\tau-1},
\end{align*}
where $D_\tau = \frac{\partial \mathbf{h}_{\tau}}{\partial \mathbf{h}_{\tau-1}}$, and $G_\tau = \frac{\partial \mathbf{h}_{\tau}}{\partial \boldsymbol{\phi}_{\tau-1}}$.
In this way, the complexity of learning no longer depends on the temporal horizon. However, this comes at the cost of memory requirements $\mathcal{O}(|\boldsymbol{\phi}||\mathbf{h_t}|)$, which makes RTRL rarely used in ML practice.

\section{Experiments}
We conduct empirical studies to assess the proposed algorithm's performance and computational cost. We compare them to widely used hyperparameter optimization methods, on standard learning problems such as MNIST and CIFAR10 classification.  To better understand the properties of OHO, we examine hyperparameter dynamics during training and how they are affected by various meta-hyperparameter settings. We explore layerwise hyperparameter sharing as a potential way to trade-off between computational costs and the flexibility of the learning dynamics.
Lastly, we analyze the stability of the optimization procedure with respect to meta-hyperparameters.

For the MNIST~\citep{LeCun1998MNIST} experiments, we use 10,000 out of the 60,000 training examples as the validation set, for evaluating OHO outer-loop gradients. We assess two architectures: a 4-layer neural network and a 4-layer convolutional neural network, with ReLU activations and $5\times5$ kernel size. We use 128 units and 128 kernels at each layer for both networks.
For the CIFAR10 dataset~\citep{Krizhevsky2009CIFAR}, we split the images into 45,000 training, 5,000 validation, and 10,000 test images and normalize
them such that each pixel value ranges between [0, 1].  We use the ResNet18 architecture~\citep{He2016ResNet} and apply random cropping and random horizontal flipping for data augmentation \citep{Shorten2019}. For both MNIST and CIFAR10,  the meta optimizer is SGD, with meta learning rate $0.000005$ and initial weight decay coefficient $0$. We set initial learning rates to $0.001$ and $0.01$, and validation batch size to $100$ and $1000$ for MNIST and CIFAR10, respectively.

\subsection{Performance}
We compare OHO against hyperparameter optimization by uniform random search and Bayesian optimization, in terms of
 the number of training runs and the computation time required to achieve a predefined level of generalization performance (test loss $\leq0.3$, see Fig.~\ref{fig:performance}).
 We vary the optimizer $g$ (Eq.~\ref{eq:update_rule}) by considering five commonly used learning rate schedulers: 1) SGD with fixed learning rate (`Fixed'), 2) Adam, 3) SGD with step-wise learning rate annealing (`Step'), 4) exponential decay (`Exp'), and 5) cosine (`Cosine') schedulers \citep{Loshchilov2017}.

 We optimize the learning rate and weight decay coefficient, and all models are trained for 100 and 300 epoch for MNIST and CIFAR10, respectively.
Additionally, we tune step size, decay rate, and momentum coefficients for Step, Exp, and Adam, respectively.
In our experiments, we find that OHO usually takes a single run to find a good solution.
The single training run takes approximately $12\times$ longer than training the network once without OHO, and yet, the other hyperparameter algorithms require multiple training of neural network to reach the same level of performance.
Overall OHO ends up significantly faster at finding a good solution in terms of total wall-clock time (Fig.~\ref{fig:performance}, bottom row).

For each method, we measure the final test loss distribution obtained for various initial learning rates and weight decay coefficients.
The initial learning rates and weight decay coefficients were randomly chosen from the range of $[0.0001,0.2]$ and $[0,0.0001]$, respectively.
For both MNIST and CIFAR10, our procedure, which dynamically tunes the learning rate and degree of regularization during model learning, results in better generalization performance and much smaller variance compared to other optimization methods (Fig.~\ref{fig:perm_stability}, `global OHO').
This demonstrates the robustness of OHO to the initialization of hyperparameters and makes it unnecessary to train multiple models in order to find good hyperparameters.
These results also suggest that there may be a systematic advantage in jointly optimizing parameters and hyperparameters.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth,page=3]{figs/figsOHO.pdf}
    \vspace{-0.5cm}
    \caption{(A) The learning rate dynamics for each layer during the early stage of the training for different initializations.
    (B) Test loss comparison for layerwise vs.\ global OHO, for several initial learning rates.}
    \label{fig:layewise_lr}
    \vspace{-0.25cm}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\linewidth,page=4]{figs/figsOHO.pdf}
    \vspace{-0.25cm}
    \caption{Test loss and wallclock time statistics when using OHO to optimize different numbers of hyper-parameter sets, where each set contains three hyper-parameters:
    two learning rates for weights and bias, and one L2 weight decay coefficient. The hyper-parameter sets are allocated to neighboring layers to evenly partition the full network. Layerwise OHO has 6 sets for MNIST, 62 for CIFAR10. }
    \label{fig:tradeoff_perm_speed}
    \vspace{-0.25cm}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics [width=0.9\linewidth,page=7]{figs/figsOHO.pdf}
    \vspace{-0.25cm}
    \caption{
    (A) Learning rate dynamics, after being manually re-initialized to 0.2 at epoch 30, 50, 90, 180, and 270 (vertical dashed lines).
    (B) Test loss of dynamic learning rate (solid line) and non-corrected learning rate (dashed line) after the reset.
    (C) Weight decay dynamics, after being manually re-initialized, details as in (a).
    (D) Corresponding test loss on CIFAR10.}
    \label{fig:perm_corruption}
    \vspace{-0.5cm}
\end{figure}

\subsection{Layerwise OHO}
In the earlier experiment (Fig.~\ref{fig:perm_stability}, `layerwise OHO'), we saw that {\em layer-specific} learning rates and weight decay coefficients can lead to even better test loss on CIFAR10. Here we further study richer hyperparametrizations.
To do so, we examine the effects of layerwise OHO on final test performance and quantify the trade-offs between performance and computational time, as we vary the number of hyperparameters.
First, we experiment with having a separate learning rate and weight decay coefficient per layer. We train both global OHO and layerwise OHO on CIFAR10 with several initial learning rates, $\alpha_0 = [0.1,0.15,0.2]$, and no initial regularization, $\lambda_0 = 0$, and compare their test losses (Fig.~\ref{fig:layewise_lr}). We find that the layerwise OHO test loss (green curves) is generally lower than that of global OHO (purple curves) is. Moreover, when analyzing the early stage of training across runs with different initial learning rates, we obverse that learning rates tend to cluster according to their layers regardless of the starting points, and their values are layer-specific (Fig.~\ref{fig:layewise_lr}).  This suggests that layerwise grouping of the hyperparameters may be an appropriate thing to do in such models.

We analyze the performance and computational speed trade-off of layerwise OHO with coarser hyperparmater grouping schemes.
If we define one hyperparameter set  as the learning rates for the weight and bias, together with the L2 regularization weights, then global OHO has one set, while layerwise OHO has as many sets as there are layers in the network (6 for MNIST and 62 for CIFAR10). We can then interpolate between these extremes. Specifically, for $k$ hyperparameter sets, we partition the network layers into $k$  neighboring groups, each with their own hyperparameter set (Fig.~\ref{fig:tradeoff_perm_speed}). For example, for $k=2$, the 6-layer network trained on MNIST would be partitioned into top 3 and bottom 3 layers, with 2 separate hyperparameter sets. For each configuration, we ran OHO optimization 10 times with different random initializations.
As shown in Fig.~\ref{fig:tradeoff_perm_speed}, the average performance improves as the number of hyperparameter sets increases, while the variance shrinks. In contrast, the wallclock time increases with the number of hyperparameter sets. This demonstrates a clear trade-off between the flexibility of the learning dynamics (and subsequently generalization performance) and computational cost, which argues for selecting the richest hyperparameter grouping scheme possible given a certain computational budget.

\subsection{Response to hyperparameter perturbations}

Due to its online nature, the OHO algorithm lends itself to being used for non-stationary problems, for instance when input statistics or task objectives change over time. Here we emulate such a scenario by perturbing the hyperparameters directly under a static learning scenario, and track the evolution of (meta-)learning in response to these perturbations.
In particular, we initialize the learning rates to a fixed value (0.1) and then reset them to a higher value (0.2) at various points during CIFAR10 training (Fig.~\ref{fig:perm_corruption}A).
We observe that the learning rates decrease rapidly early on and that after being reset, they drop almost immediately back to the value before the perturbation. This illustrates  the resilience and rapid response of OHO against learning rate corruption which reflects sudden change in an environment. We further compare the performance of OHO to the setting where the learning rate is fixed after re-initialization. Figure~\ref{fig:perm_corruption}B shows the test losses for fixed (dashed line) and OHO (solid line). We find that the dynamic optimization of hyperparameters leads to systematically better performance.

In the next set of experiments we reset the weight decay coefficient to zero at various epochs: 30, 50, 90, 180, and 270.
The results are similar, with  the weight decay coefficients quickly adapting back to their previous values (Fig.~\ref{fig:perm_corruption}C and D). The test loss fluctuates less for continually optimized OHO, relative to the fixed hyperparameters setting.
Altogether, these results suggest that joint optimization of parameters and hyperparameters is more robust to fluctuations in a learning environment.

\begin{figure}
        \centering
    \includegraphics[width=\linewidth,page=6]{figs/figsOHO.pdf}
    \vspace{-0.5cm}
    \caption{Sensitivity analysis on CIFAR10. A) The comparison of learning when the outer-gradient $\frac{\partial L_{\mathcal{D}}}{\partial \boldsymbol \varphi}$ is subject to training vs validation dataset. B) The performance when using 100, 500, and 1000 validation batch sizes. }
    \label{fig:sensitivity_analysis}
    \vspace{-0.25cm}
\end{figure}

\begin{SCfigure}
\includegraphics[width=0.5\linewidth,page=10]{figs/figsOHO.pdf}
\vspace{-0.5cm}
     \caption{Test loss as a function of learning time, when the influence get reset to zero at every 1, 100, and 1000 steps; `no reset' corresponds to standard infinite-horizon global OHO. }
     \label{fig:sensitivity_analysisC}
\end{SCfigure}

\subsection{Long-term dependencies in learning and their effects on hyperparameters }
Past work has suggested that the temporal horizon for the outer-loop learning influences both performance and the stability of hyperparameter optimization \citep{Metz2019}. In our case, the influence matrix is computed recursively over the entire learning episode and there is no strict notion of a horizon. We can however control the extent of temporal dependencies by resetting the influence matrix at a predefined interval (effectively forgetting all previous experience). We compare such resetting against our regular infinite-horizon procedure (global-OHO; see Fig.~\ref{fig:sensitivity_analysisC}). We find that both the average test loss and its variance decrease as we reset less frequently; ultimately global OHO performs the best.
Thus, we conclude it is important to take into account past learning dynamics for improving generalization performance.

\subsection{Sensitivity analysis for meta-hyperparameters.}

In order to better understand how OHO works in practice, we explore its sensitivity to meta-hyperparameter choices: the validation dataset, initial learning rates, and meta learning rates.

\paragraph{Validation dataset.}

We explore the gradient of the meta-optimization objective $\nabla_{\boldsymbol \varphi} \mathcal{L}_{\mathcal{D}_{\text{val}}}$ in Eq.~\ref{eq:oho_gradient}, where the outer gradient $\frac{\partial \mathcal{L}_{\mathcal{D}_{\text{val}}}}{\partial {\boldsymbol \theta}_{\tau+1}}$ is computed with respect to the validation dataset $\mathcal{D}_{\text{val}}$. To show the importance of having a separate validation loss for the outer loop optimization, we replace the validation dataset with the training dataset, and find that it yields to systematically worse generalization performance (Fig.~\ref{fig:sensitivity_analysis}B).
Expectedly, hyperparameters overfit to the training loss when using the training dataset, but not when using the validation dataset. When using the validation dataset, OHO prevents overfitting by shrinking the learning rates and increasing the weight regularization over time. In other words, OHO early-stops on its own.

In some practical applications it may prove computationally cumbersome to use the full validation dataset. We experiment using a random subset of validation set to compute the outer gradient while varying the validation minibatch size: 100, 500, and 1,000. Unsurprisingly we find that bigger is better (Fig.~\ref{fig:sensitivity_analysis}C): the test loss is lowest for the largest size considered, while the test loss for the smallest (100) minibatch size is on par with the step-wise learning rate scheduler. The performance gradually changes as we vary the size, which allows us to make trade-off between computational cost and the final performance.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth,page=9]{figs/figsoho.pdf}

    \vspace{-0.5cm}
            \caption{Norm of influence matrix for different meta learning rates and initial learning rates. (A) Large meta learning rates lead to a discontinuous landscape (red curve). (B) The initial learning rates $1e^{-3}$ and $1e^{-2}$ lead to instability when combined with large meta learning rates (green and blue curve).}

    \label{fig:normJ_fix_mlr}
    \vspace{-0.5cm}
\end{figure}

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.75\linewidth,page=8]{figs/figsOHO.pdf}
    \vspace{-0.5cm}
    \caption{(A,B) Demonstrating that the influence matrix norm w.r.t the learning rate and the weight coefficient plateaus, rather than exploding. (C) The gradients within a single epoch (100 updates) do not correlate with each other. The average and the standard deviation of the correlation decrease over time.}
    \label{fig:stability}
    \vspace{-0.25cm}
\end{figure*}

\iffalse
\begin{figure*}[t]
    \centering
    \begin{minipage}{0.24\textwidth}
    \includegraphics[width=\linewidth]{figs/sequential_mlp_tr.pdf}
    \vspace{-0.5cm}
    \subcaption{Training Loss}
    \label{fig:cont_mlp_tr}
    \end{minipage}
    \begin{minipage}{0.24\textwidth}
    \includegraphics[width=\linewidth]{figs/sequential_mlp_te.pdf}
    \vspace{-0.5cm}
    \subcaption{Test Loss}
    \end{minipage}
    \begin{minipage}{0.24\textwidth}
    \includegraphics[width=\linewidth]{figs/sequential_mlp_lr.pdf}
    \vspace{-0.5cm}
    \subcaption{Learning Rate}
    \end{minipage}
    \begin{minipage}{0.24\textwidth}
    \includegraphics[width=\linewidth]{figs/sequential_mlp_l2.pdf}
    \vspace{-0.5cm}
    \subcaption{Weight Decay Coefficient}
    \end{minipage}
    \caption{The performance on continual learning for various optimization methods.
        We apply multiple augmentation over the MNIST dataset during 100 epochs of training.
        First 10 epoch of regular training and then 5 epochs of horizontal flips, followed by
         5 epochs of regular training and then 15 epochs of clockwise rotation,
        10 epochs of regular training and then 5 epochs of vertical flip,
         5 epochs of regular training and then 10 epochs of target class shift by one,
        10 epochs of regular training and then 10 epochs of clockwise rotation, and
        finally 10 epochs of regular training
        (10 regular, 5 hflip, 5 regular, 15 rotation, 5 vflip, 5 regular, 10 target shift, 10 regular0, 10 rotation, 10 regular).}
    \label{fig:cont_learning}
\end{figure*}

\begin{figure*}[t]
    \centering
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/dPdLrs_all.pdf}
    \vspace{-0.5cm}
    \subcaption{Influence matrix w.r.t $\alpha$}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=\linewidth]{figs/dPdL2s_all.pdf}
    \vspace{-0.5cm}
    \subcaption{Influence matrix w.r.t $\lambda$}
    \end{minipage}
    \begin{minipage}{0.32\textwidth}
    \includegraphics[width=0.68\linewidth]{figs/gradient_correlation_sgd.pdf}
    \vspace{-0.5cm}
    \subcaption{Gradients correlations}
    \end{minipage}
    \caption{Stability of the influence matrix norm for learning rates (A) and weight regularizer (B) for multiple runs with different initialization, which corresponds to multiple lines. (C) The average and the standard deviation of the norm correlations over learning.}
    \label{fig:stability}
\end{figure*}
\begin{figure*}[ht]
    \centering
    \begin{minipage}{0.48\textwidth}
    \begin{minipage}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_bag_dPdLlr_wrt_lr_curve_comp_sgd.pdf}

    \end{minipage}
    \begin{minipage}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figs/mlp_bag_dPdLl2_wrt_l2_curve_comp_sgd.pdf}

    \end{minipage}
    \subcaption{sensitivity w.r.t initial learning rate and weight coefficient}
    \label{fig:normJ_fix_mlr}
    \end{minipage}
    \begin{minipage}{0.48\textwidth}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_bag_dPdLlr_wrt_mlr_curve_comp_sgd.pdf}
        \subcaption*{Learning rates}
        \end{minipage}
        \begin{minipage}{0.47\textwidth}
        \includegraphics[width=\linewidth]{figs/mlp_bag_dPdLl2_wrt_mlr_curve_comp_sgd.pdf}
        \subcaption*{L2 Weight coefficients}
        \end{minipage}
        \subcaption{sensitivity w.r.t meta learning rate}
    \end{minipage}
    \vspace{-0.25cm}

        \caption{Norm of influence matrix w.r.t learning rates. (A) The large meta learning rate results in discontinuous norm surfaces (red curve). (B) The initial learning rates 1e-3 and 1e-2 are sensitivity when combined with large meta learning rates (green and blue curve). }
    \label{fig:stability}
    \vspace{-0.25cm}
\end{figure*}
\fi

\paragraph{Meta learning rates and initial hyper-parameters.}

\label{exp:gradient_stability}

To deploy OHO in practice, one needs to choose a meta learning rate and the initial hyperparameters. We investigate the sensitivity of OHO to these choices in order to inform users on what may be a sensible range for each meta-hyperparameter. We define the sensible range to be the region where training is stable. We consider learning to be stable when the gradients are well-defined along the entire learning trajectory, in which computing the gradient hinges upon the product of Hessians  induced by the recursion in the influence matrix. In sum, sensible meta-hyperparameters should keep the norm of the influence matrix bounded throughout the training.

We visualize the evolution of the norms of the influence matrices for the learning rates, $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\alpha}\|^2_F$, and weight decay coefficients, $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\lambda}\|^2_F$, for different meta learning rates and initial learning rates (Fig.~\ref{fig:normJ_fix_mlr}A, B). The norm is numerical ill-conditioned (blows up) with large meta learning rates ($>1e^{-4}$) at the end stage of training, but is smooth with smaller meta learning rates.
Recall that the OHO-optimized learning rates decrease during training to avoid overfitting. Later in training this eventually leads to a regime where the meta learning rate is larger than the learning rate. Some elements of the scaled gradient can then become greater than the corresponding learning rate $\alpha_{t-1}^i < \epsilon \left[\nabla_\alpha L({\boldsymbol \theta})\right]^i$, causing instability. The ideal meta learning rate should thus be small.

We further investigate the norm of influence matrix during  training. We train multiple models with different initialization of hyperparameters while fixing the meta learning rate to $5 \times 10^{-5}$.  We find that the norm rarely explodes in practice (Fig.~\ref{fig:stability} A, B). The norms plateau as training converges, but they do not explode.
We observe similar results for the norm of the Hessian matrix $\|\frac{\partial {\boldsymbol \theta}^{(T)}}{\partial {\boldsymbol \theta}^{(T-1)}}\|^2_2$ as well.

In order to understand why the norm of influence matrix is stable along the trajectory, we look into exploding gradient phenomena.
As we discussed in Section~\ref{sec:gradient_computation}, exploding gradients likely happen when a series of gradients are correlated. We compute the moving average and standard deviation of gradient correlations with the window size of 100 updates (a single epoch).
The mean correlation very rapidly approaches zero, with the standard deviation also decreasing as training progresses (Fig.~\ref{fig:stability}C). The gradients only correlate at the beginning of  training but quickly decorrelate as learning continues. This causes the norm of the influence matrix to plateau and prevents the gradients from exploding.

\iffalse
Let us analyze the $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\alpha}\|^2_2$ and $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\lambda}\|^2_2$ landscape to understand the stability of optimization.
First, we examine the norm of influence matrix over the range of learning rates $\alpha$ and weight decay coefficient $\lambda$ at every 10 epoch out of 100 epochs conditioned on three different meta learning rates, $5e^{-5}$, $1e^{-5}$, and $5e^{-6}$ (see Figure~\ref{fig:normJ_fix_mlr}A).
The curves are widely spread in the early stage of training (each curve corresponds to every 10th epoch), and then converges (highly concentrated) as training continues.
We observe that the surface blows up when the meta learning rate is $5e^{-5}$ for both $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\alpha}\|^2_2$ and $\|\frac{d{\boldsymbol \theta}^{(T)}}{d\lambda}\|^2_2$ near the final stage of training. However, the surfaces are smooth and steady for smaller meta learning rates, $1e^{-5}$ and $5e^{-6}$.

From Figure~\ref{fig:normJ_fix_mlr}b, the landscape blew up for learning rates $1e^{-3}$ and $1e^{-2}$. This time, we examine the norm of influence with respect to the range of meta learning rates given that the networks are trained using learning rates $1e^{-3}$ and $1e^{-2}$ in Figure~\ref{fig:normJ_fix_mlr}. The norm surface is less stable as the meta learning rates become larger. This illustrates that the gradient computation, which is subject to the influence matrix, depends on meta learning rate size. After seeing the results, it is not surprising that the optimization is unstable for large meta learning rate at near finally stage of training. The learning rate decreases over the training and actually becomes zero by the end stage when trained with OHO. This causes the meta learning rate to become larger than the learning rate at some point. Then, the scaled gradient can potentially  become larger than the learning rate $\alpha_{t-1} < \epsilon \nabla_\alpha \frac{d l({\boldsymbol \theta}^{(T)})}{d\alpha}$, which causes instability in training. Therefore, we want the ideal meta learning rate to be small.
\fi

\section{Discussion}
Truly online hyperparameter optimization remains an open problem in machine learning. In this paper, we presented a novel hyperparameter optimization algorithm, OHO, which takes advantage of online temporal credit assignment in RNNs to jointly optimize parameters and hyperparameters based on minibatches of training and validation data. This procedure leads to robust learning, with better generalization performance than competing offline hyperparameter optimization procedures. It is also competitive in terms of total wallclock time.

The dynamic interaction between parameter and hyperparameter optimization was found to not only improve test performance but to also reduce variability across runs. Beyond automatic shrinking of learning rates that avoids overfitting, OHO
quickly adapts the hyperparameters to compensate for sudden changes, such as perturbations of the hyperparameters and allows the learning process to reliably find better models.

The online nature of OHO updates makes it widely applicable to both stationary and nonstationary learning problems.  We expect the same set of benefits to apply to problems such as life-long learning \citep{German2019,Kurle2020}, or policy learning in deep reinforcement learning \cite{Padakandla2019,Xie2020,Igl2020}, in which the statistics of data and learning objectives change over time. Up to now, it has been virtually impossible to do automated hyperparameter tuning in such challenging learning settings. OHO on the other hand may prove to be a key steppingstone for achieving robust automated hyperparameter optimization in these domains.

In machine learning, an enormous amount of time and energy is spent on hyperparameter tuning. Out of the box, RTRL-based OHO is already quite efficient, with memory requirement that is linear in the number of parameters and outer-loop gradients computation having a similar cost to that of computing inner gradients. This covers many practical use cases of meta-optimization. We believe that OHO and related ideas can dramatically reduce the burden of tuning hyperparameters and becoming a core component of future off-the-shelf deep learning packages.

\section{Online hyperparameter optimization}
\label{sec:oho}

With our analogies in place, we are ready to construct our online hyperparameter optimization (OHO) algorithm.
We adapt the hyperparameters at each time step, such that both the parameters and hyperparameters are jointly optimized in a single training process (Fig.~\ref{fig:my_label}).

To achieve this, we update the parameters using training data and update the hyperparameters using validation data.
Then, the update rules for $\boldsymbol{\theta}$ and $\boldsymbol{\varphi}$ are:
\begin{align}
    \boldsymbol{\theta}_{\tau+1} &= \boldsymbol{\theta}_{\tau} - \alpha_{\tau} \Delta_{\mathcal{D}_{\text{tr}}} (\boldsymbol{\theta}_{\tau}) + \alpha_{\tau}  w(\boldsymbol{\theta}_{\tau},\lambda_{\tau})
    \label{eq:param_updates}
    \\
    \boldsymbol{\varphi}_{\tau+1} &= \boldsymbol{\varphi}_{\tau} - \eta \Delta_{ \mathcal{D}_{\text{val}}} (\boldsymbol{\varphi}_{\tau}),
    \label{eq:meta_update}
\end{align}
where $\Delta_{ \mathcal{D}_{\text{tr/val}}}$ is a descent step with respect to training or validation data, respectively, and $w(\boldsymbol{\theta}_{\tau};\lambda_\tau)$ is a regularization function.
We can use any differentiable stochastic optimizer, such as RMSProp and ADAM, to compute the descent direction $\Delta$.
Without the loss of generality, from now on, we use SGD to compute the $\Delta$ and use weight decay penalty as a regularizer for the rest of the paper. Expanding $\Delta(\boldsymbol{\varphi}_{\tau}; \mathcal{D}_{\text{val}})$ in Eq.~\ref{eq:meta_update}, we have
\begin{align}
     &\Delta_{\mathcal{D}_{\text{val}}} (\boldsymbol{\varphi}_{\tau}) = \nabla_{\boldsymbol{\varphi}} \mathcal{L}_{\mathcal{D}_{\text{val}}}(\boldsymbol{\theta}_{\tau+1}, \boldsymbol{\varphi}_{\tau}) \nonumber \\
                &= \nabla_{\boldsymbol{\varphi}} \mathcal{L}_{\mathcal{D}_{\text{val}}} \big(\boldsymbol{\theta}_{\tau} -\alpha_{\tau} \big( \Delta_{\mathcal{D}_{\text{tr}}}(\boldsymbol{\theta}_{\tau})- w(\boldsymbol{\theta}_{\tau},\lambda_{\tau})\big), \boldsymbol{\varphi}_{\tau}\big) \label{eq:delta_metaopt}.
\end{align}

The gradient of hyperparameters depends on $\varphi_0\cdots \varphi_\tau$ at iteration $\tau$.
We apply RTRL to compute this gradient in an online fashion.
Let us re-write the gradient expression of RTRL in Eq.~\ref{eq:rtrl} taking into account our mappings,
\begin{align}
     \nabla_{\boldsymbol{\varphi}} \mathcal{L}_{\mathcal{D}_{\text{val}}} &= \sum^T_{\tau=1}  \frac{\partial \mathcal{L}_{\mathcal{D}_{\text{val}}}}{\partial \boldsymbol{\theta}_{\tau+1}} \Bigg( \underbrace{ \frac{\partial \boldsymbol{\theta}_{\tau+1}}{\partial \boldsymbol{\theta}_{\tau}}\frac{\partial \boldsymbol{\theta}_{\tau}}{\partial \boldsymbol{\varphi}} +\frac{\partial \boldsymbol{\theta}_{\tau+1}}{\partial \boldsymbol{\varphi}_{\tau}} }_{\Gamma_{\tau+1}} \Bigg).

     \label{eq:oho_gradient}
\end{align}

Generally, meta-optimization involves computing the gradient through a gradient, as shown in Eq.~\ref{eq:delta_metaopt}.
This causes the temporal derivative to contain the Hessian matrix,
\begin{align*}
    \frac{\partial \boldsymbol{\theta}_{\tau+1}}{\partial \boldsymbol{\theta}_{\tau}} = I - \alpha_\tau H_\tau - 2 \alpha_\tau\lambda_\tau,

\end{align*}
where $H_\tau=\mathbb{E}_{\mathcal{B}}[\nabla^2_{\boldsymbol{\theta}} \mathcal{L}]$ is the Hessian of the minibatch loss with respect to the parameters.
Then, we plug $\frac{\partial \boldsymbol{\theta}_{\tau+1}}{\partial \boldsymbol{\theta}_{\tau}}$ into the influence matrix's recursive formula,
\begin{align}

    \Gamma_{\tau+1} = \big(I-\alpha_\tau H_{\tau} -2\alpha_\tau\lambda_\tau \big) \Gamma_\tau + G_{\tau}.
    \label{eq:influence_formula}
\end{align}
By approximating the Hessian-vector product
using the finite difference method,\footnote{
See \url{https://justindomke.wordpress.com/2009/01/17/hessian-vector-products/} for details.
}
we compute the gradient $\frac{d l(\boldsymbol{\theta}_{\tau+1})}{d\alpha_\tau}$ in linear time.
The formulation shows that the influence matrix is composed of all the relevant parts of the learning history and is updated online as new data come in.
In other words, it accumulates the long-term influence of all the previous applications of the hyperparameters $\boldsymbol{\varphi}$ on the parameters $\boldsymbol{\theta}$.

Overall RTRL is more efficient than truncated BPTT for hyperparameter optimization, in both computation and memory complexity.
First, RTRL greatly reduces the computational cost since the influence matrix formulation circumvents the explicit unrolling of temporal dependencies. Our approach is intrinsically online, while BPTT is strictly an offline algorithm (although it can sometimes be loosely approximated to mimic online methods, see~\citep{Lorraine2018}). Moreover, the standard memory challenge of training RNNs with RTRL does not apply in the case of hyperparameter optimization, as the number of hyperparameters is in general much smaller than the number of parameters. For example, one popular choice of $\varphi$ is to just include a global learning rate and some regularization coefficients.

\subsection{Gradient computation}
\label{sec:gradient_computation}
Since we frame hyperparameter optimization as RNN training, one might reasonably wonder about the fundamental issues of vanishing and exploding gradients~\citep{hochreiter2001gradient} in this setup. Gradients do not vanish, because the SGD update in Eq.~\ref{eq:sgd} is additive. This resembles RNN architectures such as long- short-term memory~\citep[LSTM;][]{hochreiter1997long} and gated recurrent units~\citep[GRU;][]{chung2014empirical}, both of which were explicitly designed to have an additive gradients, as a way of avoiding the vanishing gradient issue. Nonetheless, exploding gradient may still be an issue and is worth further discussion.

\begin{figure}
    \centering
    \includegraphics[page=5,width=0.75\linewidth]{figs/figsOHO.pdf}
    \vspace{-0.25cm}
    \caption{Overview of OHO joint  optimization procedure. The parameters and hyperparameters update in parallel based on training/loss gradients computed using the current minibatch.}
    \label{fig:my_label}
    \vspace{-0.5cm}
\end{figure}

Because the influence matrix $\Gamma_{\tau+1}$ is defined recursively, the gradient of the validation loss with respect to ${\boldsymbol \varphi}$ contains the product of Hessians:
\begin{align*}
{\small
    \nabla_{\boldsymbol{\varphi}} \mathcal{L}_{\mathcal{D}_{\text{val}}} = \Bigg\langle G_{\tau}, -\sum^{\tau}_{i=0} \Bigg(\prod^{\tau}_{j=i+1} (I-\alpha_j H_j-2\alpha_j\lambda_j)\Bigg) G_i\Bigg\rangle ,

    }
\end{align*}
where $\langle \cdot, \cdot \rangle$ denotes an inner product. $G_i$ and $H_j$ are the gradient and Hessian of the loss $\mathcal{L}(\boldsymbol{\theta}_\tau)$ at iteration $i$ and $j$, respectively. The product of the Hessians can lead to gradient explosion, especially when consecutive gradients are correlated~\citep{Metz2019}.
It remains unclear whether the gradient explosion is an issue for RTRL (or, more generally, forward-differentiation based optimization). We study this in our experiments.

\iffalse
In order to analyze the exploding gradient issue, we consider the gradient of the validation loss after $\delta$ steps w.r.t. the learning rate at time $t$:
\begin{align*}
    \frac{\partial L(\boldsymbol{\theta}_t+\delta)}{\partial \eta^{(t)}}
    &=
    \frac{\partial L(\boldsymbol{\theta}_t+\delta)}{\partial \boldsymbol{\theta}_{t+\delta}}
    \frac{\partial \boldsymbol{\theta}_{t+\delta}}{\partial \boldsymbol{\theta}_t}
    \frac{\partial \boldsymbol{\theta}_t}{\partial \eta^{(t)}}

\end{align*}
where $\eta^{(t)}$ is the use of the learning rate $\eta$ at time $t$, and $g(\theta) = \frac{\eta}{B} \sum_{(x,y) \in B} \nabla_\theta l(y, f(x; \theta))$. More specifically, we look at the temporal derivative:
\begin{align*}
    \frac{\partial \theta_{t+\delta}}{\partial \theta_t}
    =
    \frac{\partial \theta_{t+\delta}}{\partial \theta_{t+\delta-1}}
    \cdots
    \frac{\partial \theta_{t+1}}{\partial \theta_{t}},
\end{align*}
where the one-step temporal derivative is written as
\begin{align}
\label{eq:1-step-temporal-derivative}
    \frac{\partial \theta_{t+1}}{\partial \theta_{t}}
    =
    I - \eta^{(t+1)} H(\theta_t).
\end{align}
We use $H(\theta_t)$ as a shorthand for the Hessian of the mini-batch loss w.r.t. the parameters:
\begin{align}
\label{eq:hessian}
    H(\theta_t) = \frac{1}{B} \sum_{(x,y) \in B} \text{Hess}_{\theta} l(y, f(x; \theta_t)).
\end{align}
\fi

\section{Introduction}

The success of training complex machine learning models critically depends on good choices for hyperparameters that control the learning process. These hyperparameters can specify the speed of learning as well as model complexity, for instance, by setting learning rates, momentum, and weight decay coefficients. Oftentimes, finding the best hyperparameters requires not only extensive resources, but also human supervision.

Well-known procedures, such as random search and Bayesian hyperparameter optimization, require fully training many models to identify the setting that leads to the best validation loss~\citep{Bergstra2012, Snoek2012}. While popular for cases when the number of hyperparameters is small, these methods easily break down when the number of hyperparameters increases. Alternative gradient-based hyperparameter optimization methods take a two level approach (often referred to as the `inner' and the `outer loop'). The outer loop finds potential candidates for the hyperparameters, while the inner loop is used for model training. These methods also require unrolling the full learning dynamics to compute the long-term consequences that perturbation of the hyperparameters has on the parameters~\citep{luketina2016scalable, Lorraine2018, Metz2019}. This inner loop of fully training the parameters of the model for each outer loop hyperparameter update makes existing methods computationally expensive and inherently offline.
It is unclear how to adapt these ideas to nonstationary environments, as required by advanced ML applications, e.g.\ lifelong learning~\cite{German2019,Kurle2020}.
Moreover, while these type of approaches can be scaled to high-dimensional hyperparameter spaces~\citep{Lorraine2018}, they can also suffer from stability issues~\citep{Metz2019}, making them nontrivial to apply in practice.
In short, we are still missing robust and scalable procedures for online hyperparameter optimization.

We propose an alternative class of online hyperparameter optimization (OHO) algorithm that update hyperparameters in parallel to training the parameters of the model. At the core of our framework is the observation that hyperparameter optimization entails a form of \emph{temporal credit assignment}, akin to the computational challenges of training  recurrent neural networks (RNN). Drawing on this similarity allows us to adapt real-time recurrent learning (RTRL)~\citep{Williams1989}, an online alternative to backprop through time (BPTT), for the purpose of online hyperparameter optimization.

We empirically show that our joint optimization of parameters and hyperparameters yields systematically better generalization performance compared to standard methods. These improvements are accompanied sometimes with a substantial reduction in overall computational cost. We characterize the behavior of our algorithm in response to various meta-parameter choices and to dynamic changes in hyperparameters. OHO can rapidly recover from adversarial perturbations to the hyperparameters. We also find that OHO with layerwise hyperparameters provides a natural trade-off between a flexible hyperparameter configuration and computational cost. In general, our framework opens a door to a new line of research, that is, real-time hyperparameter learning.

Since we have posted the initial version of this preprint, it has been brought to our attention that the early work by~\citep{Franceschi2017} presented a very closely related method to ours, from a different perspective. \citet{Franceschi2017} studied both forward and reverse gradient-based hyperparameter optimization. They derived RTRL \citep{Williams1989} by computing the gradient of a response function from a Lagrangian perspective. In contrast, here we directly map an optimizer to a recurrent network, by taking the process of fully training a model  as propagating a sequence of minibatches through an RNN, which allows us to use RTRL for the purpose of online hyperparameter optimization. Moreover, if \citet{Franceschi2017} included a single small experiment comparing their approach to random search on phonetic recognition, here we empirically evaluate our algorithm's strengths and properties extensively, and compare them to several state of the art methods.
In doing so, we demonstrate OHO's i) effectiveness: performance w.r.t the wall clock time improves existing approaches, ii) robustness: low variance of generalization performance across parameter and hyparparameter initializations, iii) scalability: sensible trade-off between performance vs. computational cost  w.r.t.\ the number of hyperparameters, and iv) stability of hyperparameter dynamics during training for different choices of meta-hyperparameters.

\input{background}
\input{setup}
\input{hyperopt}
\input{relatedwork}
\input{experiment}

\nocite{langley00}

\section{Related work}
Optimizing hyperparameters based on a validation loss has been tackled from various stances. Here, we provide a brief overview of state-of-the-art approaches to hyperparameter optimization.

\paragraph{Random search.}

For small hyperparameter spaces, grid or manual search methods are often used to tune hyperparameters.
For a moderate number of hyperparameters, random search can find hyperparameters that are as good or better than those obtained via
 grid search, at a fraction of its computation time \citep{Bengio2000, Bergstra2012, Bergstra2011, Jamieson2015, Li2016}.

\paragraph{Bayesian optimization approaches.}
Bayesian optimization (BO) is a smarter way to search for the next hyperparameter candidate \citep{Snoek2012, Swersky2014, Snoek2015, Eriksson2019, Kandasamy2020} by explicitly keeping track of uncertainty.
BO iteratively updates the posterior distribution of the hyperparameters and assigns a score to each hyperparameter based on it. Because evaluation of hyperparameter candidates must largely be done sequentially, the overall computational benefit from BO's smarter candidate proposal is only moderate.

\paragraph{Gradient-based approaches.}
Hyperparameter optimization with approximate gradient (HOAG) is an alternative technique that iteratively updates the hyperparameters following the gradient of the validation loss.
HOAG approximates the gradient using an implicit equation with respect to the hyperparameters \cite{Domke2012, Pedregosa2016}.
This work was extended to DNNs with stochastic optimization settings by \citet{Maclaurin2015, Lorraine2020}.  In contrast to our approach, which uses the chain rule to compute the exact gradient online, this approach exploits the implicit function theorem and inverse Hessian approximations.
While the hyperparameter updates appear online in form, the method requires the network parameters to be nearly at a stable point. In other words, the network needs to be fully trained at each step.

\citet{Metz2019} attempt to overcome the difficulties of backpropagation through an unrolled optimization process by using a surrogate loss, based on variational and evolutionary strategies. Although their method addresses the exploding gradient problem, it still truncates
backpropagation \citep{Shaban2019},
which introduces a bias in the gradient estimate  of the per step loss. In contrast, our method naturally addresses both issues: It is unbiased by algorithm design, and the gradients are empirically stable, as will be shown in Experiment~\ref{exp:gradient_stability}.

\paragraph{Forward differentiation approaches.}

Our work is closely related to \citep{Franceschi2017, Baydin2018, Donini2020}, which share the goal of optimizing hyperparameters in real-time.
\citet{Franceschi2017} first introduced a forward gradient-based hyperparameter optimization scheme (RTHO), where they compute the gradient of the response function from a Lagrangian perspective. Their solutions ends up being equivalent to computing the hypergradient using RTRL \citep{Williams1989}, as in OHO. Our conceptual framework is however different, as we derive the hyperparameter gradient by directly mapping an optimizer to a recurrent network and then applying RTRL. Unlike \citet{Franceschi2017}, whose only quantitative evaluation was performed on TIMIT phonetic recognition,  we evaluate OHO extensively in a large experimental setup and rigorously analyze its properties. \citet{Donini2020} recently proposed adding a discount factor $\gamma \in [0,1]$ to the influence function, as an heuristic designed to compensate for inner loop non-stationarities by controlling the range of dependencies of the past influence function.

One potential disadvantage of our approach is that the memory complexity grows with the number of hyperparameters. There are however several approximate variants of RTRL that can be used to address this issue.
As an example, unbiased online recurrent optimization (UORO) reduces the memory complexity of RTRL from quadratic to linear \cite{Tallec2017}.
UORO uses rank-1 approximation of the influence matrix and can be extended to rank-$k$ approximation~\citep{Mujika2018}.
Although UORO gives an unbiased estimate of the gradient, it is known to have high variance, which may slow down or impair training \cite{Cooijmans2019, Marschall2020}. It remains an open question whether these more scalable alternatives to RTRL are also applicable for online hyperparameter optimization.

\section{Learning as a recurrent neural network}

Both RNN parameter learning and hyperparameter optimization need to take into account long-term consequences of (hyper-)parameter changes on a loss. This justifies considering an analogy between the recurrent form of parameter learning (Eq.~\ref{eq:update_rule}) and RNN dynamics (Eq.~\ref{eq:rnn}).

\paragraph{Mapping SGD to a recurrent network.}
To establish the analogy, consider the following mapping:
\begin{align*}
    \text{parameters } \boldsymbol{\theta}_\tau &\rightarrow \text{state } \mathbf{h}_t \\
    \text{batch } \mathbf{B}_\tau &\rightarrow \text{input } \mathbf{x}_t\\
    \text{update rule } g(\cdot) &\rightarrow \text{recurrent function } r(\cdot)\\
    \text{hyper-parameters } \boldsymbol{\varphi} &\rightarrow \text{parameters } \boldsymbol{\phi}.
\end{align*}
Parameters ${\boldsymbol \theta}_\tau$ and recurrent network state $\mathbf{h}_t$ both evolve via recurrent dynamics, $g(\cdot)$ and $r(\cdot)$, respectively.
These dynamics are each parametrized by hyper-parameters ${\boldsymbol \varphi}$ and parameters ${\boldsymbol \phi}$.
The inputs $\mathbf{x}_t$ in the RNN correspond to the minibatches of data $\mathbf{B}_\tau$ in the optimizer, and the outputs $\mathbf{o}_t$ to
$\left\{
    f(x; {\boldsymbol \theta}_\tau)
    \Big|
    x \in B_\tau
    \right\}$, respectively.
In short, we interpret a full model training process as a single forward propagation of a sequence of minibatches in an RNN.

\paragraph{Mapping the validation loss to the RNN training loss.}

The same analogy can be made for the training objectives. We associate the per-step loss of the recurrent network to the hyperparameter optimization objective:
\begin{align*}
    \mathcal{L}(\boldsymbol{\theta}_\tau, \boldsymbol{\varphi}_\tau; \mathcal{D}_\text{val}) = \frac{1}{|\mathcal{D}_{\text{val}}|} \sum_{(\mathbf{x},\mathbf{y}) \in \mathcal{D}_{\text{val}}} l\big(\mathbf{y}, f(\mathbf{x}; \boldsymbol{\theta}_\tau(\boldsymbol{\varphi}_\tau))\big),

\end{align*}
where $l$ is the per-example loss function, with $\boldsymbol{\varphi}_\tau = \lbrace \alpha, \lambda \rbrace$  the set of tune-able hyperparameters, which in our concrete example corresponds to learning rate $\alpha$ and regularizer $\lambda$.
The hyperparameter optimization objective is calculated over the validation dataset $\mathcal{D}_{\text{val}}$.
The total loss over a full training run is $L_{\text{Total}} = \sum_{\tau=1}^{\infty} \mathcal{L}(\boldsymbol{\theta}_\tau, \boldsymbol{\varphi}_\tau; \mathcal{D}_{\text{val}}),$
where there are multiple parameter updates within a run. We use $\infty$ to emphasize that the number of epochs is often not determined {\it a priori}.

Both the hyperparameter optimization objective and the update rule, Eq.~\ref{eq:sgd}, contain the per-example loss $l$.
The former updates the parameters of RNN using $\mathcal{D}_{\text{val}}$, while the latter updates the state of RNN using training batch ${\bf{B}_{\tau}}$.
This meta-optimization loss function is commonly used for hyperparameter optimization in practice~\citep{Bengio2000, Pedregosa2016}. Nevertheless, the mapping from the optimizer to RNN and the mapping from the validation loss to RNN training objective is novel and enables us to treat the whole hyperparameter optimization process as an instance of RNN training.

redpajama_set_name: RedPajamaArXiv