\section{Introduction}

\nj{Nowadays}, 2K (1920$\times$1080) videos are \nj{widely used in} digital broadband systems. From the viewpoint of image quality, there \nj{are} two major \nj{sources of impaired quality in the current video images}. The first \nj{one} is \nj{the intrinsic} limitation of image capture systems such as lens optics, sensor resolution, focusing performances and ISO noise. These defects make \nj{inherent limitations in image quality}\sam{. Therefore} almost every videos cannot show \sam{ideal 2K image quality, with maximum frequencies.}
\nj{The} second \nj{source of corruptions is due from} image processing such as compression, resolution converting and noise reduction. \nj{Most components of a broadband system have a fixed resolution and they require} unavoidable needs of image scaling. \nj{Particularly, in order to transmit sources with lower resolution made in the past via a broadband system, video upscaling is inevitable.} Also \nj{most} of contents providers
\nj{use} heavy compression
\sam{and restoring the lost information is a difficult task \nj{which requires enhanced} techniques.}

\begin{figure}
\begin{center}
   \includegraphics[width=1\linewidth]{concept2.png}
\end{center}
\vspace{-5mm}
   \caption{\textbf{\nj{Proposed} Concept.} It would be better to learn through super-resolution-processed images. \nj{The figures are} $HR_0$ (original), $HR_1$ and $HR_3$ \nj{obtained by our method}. With the same network, \nj{(c) shows better quality than} (b). }
\vspace{-2mm}
\label{fig:concept}
\end{figure}

\nj{The early image enhancement techniques were mostly} at the level of unsharp-masking, which is \nj{just} a simple amplification \nj{operation}. With the advent of high-resolution panels, super-resolution \nj{techniques have been} developed to improve the original resolution. \nj{In} recent years, machine learning \nj{technologies are widely used} and deep learning has \nj{achieved state-of-the-art} super-resolution \nj{performance}. However, \nj{as} the resolution of the image \nj{as well as} the complexity of the algorithm increases, \nj{implementation becomes harder and harder. }
\nj{Therefore, it is necessary to develop} cheaper and more effective algorithms \nj{for image enhancement and super-resolution}.

\nj{In this paper, as shown in Figure \ref{fig:concept}, we introduce} the concept of iterative learning as an approach for image enhancement. In general, if the same \nj{algorithm} is applied repeatedly, the result will be reinforced. Especially in deep learning, we can approach \nj{this idea} without \nj{additional computational cost at inference time. What we need is additional training} with \sam{the new HR, made in the previous stage. Through this recurrent learning, the image converges to a better quality, and it \nj{is} possible to obtain \nj{an enhancement} over the existing super-resolution with a smaller-sized network. \nj{By training a small network with 200K parameters three repetitions, we could get images with better quality compared to a model with over 1M parameters.} }

\nj{The contributions of this paper are as follows:}

1. By showing an effective \nj{way of repeatedly training a simple feed-forward network for super-resolution}
through recurrent learning, \nj{we improve} the efficiency of super-resolution \nj{achieving} similar results to those \nj{of complex deep networks}
without \nj{much computational} cost.

2.	By shifting the evaluation criteria of the image quality from a simple \nj{arithmetic measure} like Mean-Squared-Error \nj{(MSE)} to a user-centered measure of \nj{VIQET \cite{VIQET}} Mean-Opinion-Score \nj{(MOS)}, \nj{which can be computed in an automated way without human involvement,} we can focus on additional \nj{subtle image quality} improvement \nj{and set the number of repetitions appropriately}.

\section{Related work}
\subsection{\nj{Image enhancement}}
Along with the development of computer vision and the spread of digital broadcasting, image processing techniques have been widely used in image production and reproduction. \nj{These techniques} can be roughly divided into \nj{noise} reduction and \nj{image} enhancement. The former includes technologies such as de-noising and de-interlacing that remove artifacts that should not be present in the original image. In the latter, sharpness and contrast \nj{are} used to amplify \nj{some} components in the original image \nj{for better image quality}. These algorithms typically perform a series of operations that widen the min-max range of the original \nj{image}, such as \nj{the} Laplacian filtering. The Laplacian filtering is inexpensive and effective, but it is not possible to perform \nj{more} complex processing \nj{such as restoring} missing texture components or broken edges.

There are several image enhancement algorithms that \nj{resolve} these weaknesses. \nj{The first one is} improving the adaptiveness \nj{of the processing region} such as \nj{Adaptive Unsharp Masking \cite{826787}, which selects the processing region wisely}. The second \nj{one} is to use two or more weights in the function to enable more complex processing, which includes \sam{Bilateral Filter \cite{4476197} and NLM Filter\cite{1467423}}. There \nj{are huge amount of} related works, but with the advent of super-resolution, research on conventional image enhancement has diminished.
{This is \nj{mainly} because the improvement is due \nj{from} the amplification of \nj{some components}, and the original frequency is maintained. Ultimately, it requires an approach such as super-resolution in order to change the shape of the original. We present a direction for developing new super-resolution \nj{images} for better \nj{quality}.}

\subsection{\nj{Super-resolution}}
Before the boom of deep learning, as a better upscaler, \nj{super-resolution algorithms were} developed. \nj{These techniques have mainly been} applied to a \nj{display} system such as a television. \nj{It} is proposed as a solution to compensate for \nj{the mismatch between the input image size and the panel resolution}.
There are two \nj{main} categories of \nj{this} technology: \nj{1)} utilizing machine learning \nj{algorithms} and \nj{2) weighted image blending using the patches} from the original \nj{image}. The former is similar to the super-resolution \nj{approaches} through the current deep learning, but \nj{with simple architectures} consisting of only one layer. The latter is a self-image generation method that retrieves the most similar patches \nj{from the surrounding areas of images after upscaling} \cite{5459271}.
Conventional super-resolution \nj{methods are} well documented in \nj{\cite{5466111}}. \nj{They} can perform upscaling while keeping the sharpness of an edge. However, there are technical limitations such as \nj{the problem for} a large upscaling ratio or impossibility of \nj{restoring} complex textures.  {\nj{These require} huge computations, but it was difficult to learn enough with a human-designed shallow network. Recently, \nj{deep neural networks have} been improved to enable more powerful super-resolution processing. This is why we use deep learning as a solution.}

\subsection{Super-resolution with deep learning}
\nj{Recently, as} a branch of deep learning \nj{research}, `single image super-resolution' has \nj{emerged}. It makes down-sized \nj{low-resolution (LR)} images and \nj{trains} convolution networks \nj{so as to reproduce the original images from the LR images}.
\nj{Compared to the revolutionarily simple} SRCNN \cite{10.1007/978-3-319-10593-2_13}, \nj{the architecture has become} heavier and more complex.
VDSR \cite{Kim_2016_CVPR} introduced VGG and ResNet into \nj{the} super-resolution filed \nj{in} 2015. ESPCN \cite{Shi_2016_CVPR} \nj{and} DRCN \cite{Kim_2016_CVPR_DRCN} present various aspects of residual structure and \nj{the} view of upscaling \nj{in} 2016. SRGAN \cite{Ledig_2017_CVPR}, SRResNET \cite{SRResnet}, DRRN \cite{Tai_2017_CVPR}, EDSR \cite{Lim_2017_CVPR_Workshops}, DenseSR \cite{Tong_2017_ICCV}, MDSR \cite{Lim_2017_CVPR_Workshops} \nj{and} Memnet \cite{Tai_2017_ICCV} \nj{developed} Resnet architecture and \nj{adjusted the} GAN solution into super-resolution \nj{in} 2017. In 2018, super-resolution is now more widely used in terms of utilization.
A good example is  \cite{Yoo_2018_CVPR}. \nj{There is a research utilizing a recursive method for super-resolution \cite{Dahl_2017_ICCV}.
\nj{Different from ours,
it tries to} \sam{produce \nj{a} face image from \nj{an} 8$\times$8 tiled mosaic
through repeated application of super-resolution, \sam{\nj{which is} similar to info-GAN \cite{NIPS2016_6399} \nj{that} produces a facial image.}

\sam{Regardless of the super-resolution algorithm,}
\nj{the general operation principle is similar. Basically, the} single image super-resolution is learned by back propagation so as to \nj{reduce} the difference between the \nj{original and the generated image through the SR network from an LR image}.

There are some limitations in these works. First, in the low-resolution images used in other SR papers, the improvement of the image \nj{quality} by \nj{increasing the resolution} is remarkable. However in a high-resolution image, the improvement caused by an increase in resolution \nj{is} not obvious. Second, the advances in super-resolution have led \nj{to an increased performance sacrificing the computational complexity.}
\nj{This high computational demand makes it hard to implement super-resolution with hardware. Therefore,} it is necessary to develop a method that can achieve better performance with a more simple network. Finally, \nj{since current super-resolution algorithms are} learned \nj{only from the} original \nj{images, they target} the original and \nj{we cannot expect} image quality beyond the original. {\nj{In the following section, we} present a new approach to get \nj{images with better quality at a low cost} through recurrent learning.}

\section{\nj{Proposed method}}

\subsection{\nj{Recurrent Training Strategy (RTS)}}
The core of this paper is the recurrent \nj{training strategy of an SR network which is achieved by} reducing the \nj{size of the output image of the SR network} to its original size and \nj{retraining the SR network with this new target}. The \nj{proposed} recurrent \nj{training strategy} is \nj{briefly} shown in Figure \ref{fig:overall}.

\begin{figure}
\includegraphics[width=\linewidth]{stage.png}
\vspace{-4mm}
\caption{\nj{\textbf{Recurrent Training Strategy(RTS).} LR is a low resolution image obtained by down-scaling the original image ($HR_0$). One stage is composed of the SR training phase (Phase A) and the image enhancement phase (Phase B). By successive application of these two phases, we can obtain a better SR network and enhanced images.} \sam{Although it takes multiple times to learn, \nj{at inference time,} it is possible to obtain \nj{images with quality equivalent to the finally learned $HR_{n}$ level} at once \nj{by feed-forwarding with the learned weights}.
} }
\label{fig:overall}
\end{figure}

\nj{As can be seen in the figure, a stage of RTS} consists of two phases. In the first phase \nj{(Phase A), SR network is trained}. \sam{In this phase, the SR is learned} so that the target image is obtained from the degraded LR image, and the weight of the super-resolution network is updated. The LR image is a \nj{low quality} image which is corrupted through down-scaling, adjusting compression \nj{and so on.}
In Stage 1\sam{, the target is set to the original image, $HR_0$ and} the network is \nj{tuned such} that the LR is improved to the image quality of the original $HR_0$ level.

The second phase, Phase B, is the process of creating a new target image using the network \nj{obtained in} Phase A. \nj{This is done by inputting the original image $HR_0$ to the trained network} of Phase A \nj{and we can expect that the output image $SR_1$ has higher quality than the original image $HR_0$. By downscaling the output image $SR_1$, we can obtain a new image $HR_1$ which acts as a target in the next stage. }
Generally, super-resolution is accompanied by increasing the size of the image to increase the image resolution. So the additional down-scaling process is required because the resolution must be the same as the original for \nj{retraining the network}.

To summarize, we have to go through two steps to get a new HR \nj{image}. The former is the process of obtaining the super-resolution \nj{function $SR(\cdot)$}, and the latter is the \nj{process of obtaining enhanced image by down-scaling the output image of the SR network back to the original size}. Then, we perform this process repeatedly on a stage basis, so that the same super-resolution processing can be \nj{more effective}.
\nj{The overall procedure of Stage $n$ can be summarized as follows:
\begin{equation}
    \begin{split}
        &\text{Phase A: } SR_n = \argmin_{SR} \frac{1}{2}||HR_{n-1} - SR(LR)||^2 \\
        &\text{Phase B: } HR_n = d(SR_n(HR_0)).
    \end{split}
\end{equation}
Here, $SR(\cdot)$ and $d(\cdot)$ denote the super-resolution network and the down-scaling operation, respectively. $LR$ is an image which has lower resolution than the original $HR_0$.
}

By applying super-resolution to the original image, a \nj{higher quality image is expected to be}  obtained. Even if the obtained image is down-scaled to the original size, \nj{components from the additionally} generated resolution \nj{may remain}, this is the same \nj{as an UHD (Untra HD) channel that looks better than the existing channels on a FHD \nj{(full HD)} TV}.
If \nj{$HR_n$} has a better image quality than \nj{$HR_{n-1}$}, we can assume that the super-resolution \nj{network} can learn from \nj{this} better target, and after \nj{the learning, the resultant network will be better.}
\nj{Therefore, assuming} that we go through this additional process in a direction that gradually improves the results, we will eventually reach the maximum improvement point of each image. And it \nj{will perform better than} the \nj{current} super-resolution \nj{network trained to mimic the original image}.

\begin{figure*}
    \includegraphics[width=0.95\linewidth]{fig_2.png}
    \caption{\textbf{Used network for recurrent learning.} {This network was developed exclusively for recursive learning systems. The biggest difference from the existing SR network is \nj{that it was} designed to produce the best result at the \sam{blue layer} \nj{where the} residual is added \nj{with only 3 channels.}}
    \nj{The model} is also very light and easy to control the number of layers.}
\vspace{-2mm}
\label{fig:network}
\end{figure*}

\subsection{\nj{Reason for enhanced quality}}
If $y$ is the original image, $U$ is the up-scaler, $D$ is the down-scaler, and $R$ is the residual component, the super resolution trained by \nj{$x=D(y)$} can be expressed as follows:
\begin{equation}
    SR(x) = U(x) + R(x),
    \label{eq:srx}
\end{equation}
\begin{equation}
    SR(D(y)) = U(D(y)) + R(D(y)).
\end{equation}
Here, we can suppose $D(U(y))$ is \nj{the} same as $y$, and $U(D(y))$ \nj{becomes} \nj{a locally blurred version of} $y$:
\begin{equation}
    D(U(y)) = y, \quad
    U(D(y)) = \bar{y}
    \label{eq:udy}
\end{equation}
\nj{Combining (\ref{eq:srx}) and (\ref{eq:udy}), we get }
\begin{equation}
    R(x) = SR(x) - \bar{y}.
\end{equation}
\nj{Because} the goal of learning is to minimize \nj{the error $||y-SR(x)||$}:
\begin{equation}
    SR^* = \argmin_{SR} \left( \frac{1}{2} \left \| y - SR(x)\right \| ^{2} \right),
\label{eq_3}
\end{equation}
\nj{if} enough learning has been done, $R(x)$ \nj{becomes} similar to $y-\bar{y}$:
\begin{equation}
    R(x) \approx y - \bar{y}.
\label{eq:approx1}
\end{equation}
In the first stage, after \nj{the} super resolution is fully learned, \nj{instead of $D(y)$,} we can input $y$ to \nj{the SR network}, \nj{then it becomes}
\begin{equation}
    SR_{1}(y) = U(y) + R_{1}(y).
\end{equation}
\nj{By down-scaling this, we get the target image for stage 2 which can be expressed as}
\begin{equation}
    HR_{1}(y) = D(SR_{1}(y)) = D(U(y) + R_{1}(y)).
\end{equation}
\nj{Assuming a linear\footnote{For a nonlinear down-scaler such as nearest neighbor or bicubic interpolation, we can take Taylor series expansion.} down-scaler, as $D(U(y)) = y$, it becomes}
\begin{equation}
    HR_{1}(y) = y + D(R_{1}(y)).
\end{equation}
\nj{The distributions of zero-mean residual components of different sized images will be approximately equal, that is $p(z|z\in R_1(x)) = p(z|z\in R_1(y))$, unless there are not much high frequency components only existing in $y$ and not in $x$. By down-scaling $R_1(y)$, the distribution will be more contracted near zero and we can approximate $D(R_1(y))$ with a contracted version of $R_1(x)$ as }
\begin{equation}
    D(R_{1}(y)) \approx \alpha R_{1}(x) , \quad ( 0 \le \alpha \le 1 ).
\label{eq:Dtoalpha}
\end{equation}
Finally, with (\ref{eq:approx1}), $HR_{1}(y)$ \nj{can} be approximated \nj{as}
\begin{equation}
    HR_{1}(y) \approx y + \alpha R_{1}(x) \approx y + \alpha (y-\bar{y}).
\end{equation}
\nj{The above equation} has high similarity \nj{with} unsharp masking. In other words, $HR_{1}(y)$ has a high probability that the \nj{components with a large deviation from the mean} are amplified. And it can be expected that \nj{the new target image} $HR_{1}(y)$ can \nj{have} improved sharpness compared to $y$.

When \nj{the} super resolution \nj{network is trained with this new target, it becomes
\begin{equation}
    HR_2(y) \approx y + \alpha R_2(x) \approx y + \alpha(y-\bar{y}) + \alpha^2 (y-\bar{y}).
\end{equation}
This is obtained by replacing $y$ in (\ref{eq:approx1}) with the new target $y+\alpha(y-\bar{y})$ in computing $R_2(x)$.
As stage goes on, it becomes}
\nj{
\begin{equation}
    HR_{n}(y) \approx y + \sum_{i=1}^n \alpha^i (y-\bar{y}).
\label{final_eq}
\end{equation}
}
\nj{Since $\alpha$} is smaller than 1, \wacvadded{(Experimental results show that the average $\alpha$ for stage 1 on DIV2K valid set is 0.779 with the standard deviation of 0.098.)} \nj{higher order $\alpha$ terms} will disappear and $HR_n$ can be expected \nj{to} convergence to \nj{a certain point $\hat{y}$ not far from $y$}.
\nj{Using the SR network, the actual difference $y-\bar{y}$ is learned by a nonlinear function \wacvadded{and it makes hard to guarantee convergence, these non-linearity} has much better performance than the simple USM algorithm. As iteration goes on, the difference tend to converges to $\hat{y}-\bar{y}$ and the image quality of $\hat{y}$ becomes better than that of the original image $y$.}

\subsection{\nj{Configuration of super-resolution network}}
We set up a convolution network based on the residual layer. The use of multiple layers of residual type has the advantage of automatically adapting image filters of various sizes. In order to make use of these advantages and to easily control the number of parameters, we constructed the \nj{network shown in Fig. \ref{fig:network}}.

{This network is designed exclusively for recurrent learning, which uses \nj{a subsequent} down-scaler \nj{after the operation of super-resolution}. The biggest difference from \nj{the} regular super resolution network \nj{is that} it produces the best output of the same resolution with original, not \nj{in an} upscaled domain. Learning proceeds in the same way as the existing super resolution, but since the expected output is synthesized in 3ch at the original resolution, optimum result is produced at the corresponding step.}

It is composed \nj{of an} input module, \nj{an} output module and three \nj{64-channel residual networks using 3$\times$3 filters}
The input module \nj{consists of} 128$\times$128 RGB (3ch) input and 64ch output. After the input module, a \nj{series of residual layers, all of which have the same number of channels (64 ch) with 3$\times$3 convolutions, is applied. As shown in the figure,} the base model uses 3 layers. In the output module, the final residual \nj{layer does not make use of} ReLU because \nj{the difference between the input and the output} will have both \nj{a positive and negative} value. {Here, the output with blue residual will be next stage's HR, because the next channel extension and pixel shuffle has no ReLU and it can be canceled by the down-scaler.} The residual layer is \nj{based} on SRResNet \cite{SRResnet} and VDSR \cite{Kim_2016_CVPR}. \nj{The final residual without ReLU is} from VDSR, and \nj{the} inner residual layer \nj{shares the} concept \nj{with} SRResNet.

Finally, we need a down-scaler to convert the \nj{up-scaled} image to the same size as the original, using the Lanczos scaler. Lanczos filter has better performance than linear \nj{and bicubic filters}, and it makes $\alpha$ greater in \nj{(\ref{eq:Dtoalpha})}.

\begin{table}
\caption{\textbf{Image difference ratio of validation set.} \nj{In our model with the proposed} recurrent \nj{training strategy}, the change of difference ratio \nj{from the previous stage is minimum at stage 4.
}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
	Stage N&Difference Ratio &Delta from (N-1)\\
	& form the Original(\%)& \\
\hline\hline
	\texttt{HR1}&0.95&0.95\\
	\texttt{HR2}&1.35&0.40\\
	\texttt{HR3}&1.47&0.12\\
	\texttt{HR4}&1.48&\textbf{0.01}\\
	\texttt{HR5}&1.51&0.03\\
	\texttt{HR6}&1.55&0.04\\
	\texttt{HR7}&1.62&0.07\\
	\texttt{HR8}&1.76&0.14\\
\hline
\end{tabular}
\end{center}
\vspace{-2mm}
\label{table:image_diff}
\end{table}

\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{fig_3.png}
\caption{\textbf{Change of image difference \nj{vs.} stage.} It can be seen that the difference converges to almost zero in HR4 (Yellow Line). But with repeated application of unsharp masking, we can observe the delta \nj{increases} linearly (Green Line). }
\label{fig:image_diff_graph}
\end{figure}

\begin{figure*}
\begin{center}
    \begin{tabular}{c}
    \hspace{-5mm}(Original)\hspace{23mm}(HR1)\hspace{27mm}(HR2)\hspace{27mm}(HR3)\hspace{27mm}(HR4)\\
    \vspace{-4mm}
    \\
    \hspace{-3mm}
    \includegraphics[width=1.35in]{0071033hr0x4resized.png}
    \includegraphics[width=1.35in]{0071033hr1x4resized.png}
    \includegraphics[width=1.35in]{0071033hr2x4resized.png}
    \includegraphics[width=1.35in]{0071033hr3x4resized.png}
    \includegraphics[width=1.35in]{0071033hr4x4resized.png}\\
    \vspace{-4mm}
    \\
    \hspace{-3mm}
    \includegraphics[width=1.35in]{0576091hr0x4resized.png}
    \includegraphics[width=1.35in]{0576091hr1x4resized.png}
    \includegraphics[width=1.35in]{0576091hr2x4resized.png}
    \includegraphics[width=1.35in]{0576091hr3x4resized.png}
    \includegraphics[width=1.35in]{0576091hr4x4resized.png}\\
    \vspace{-4mm}
    \\
    \hspace{-3mm}
    \includegraphics[width=1.35in]{0018011hr0x4resized.png}
    \includegraphics[width=1.35in]{0018011hr1x4resized.png}
    \includegraphics[width=1.35in]{0018011hr2x4resized.png}
    \includegraphics[width=1.35in]{0018011hr3x4resized.png}
    \includegraphics[width=1.35in]{0018011hr4x4resized.png}\\
    \vspace{-4mm}
    \\
    \hspace{-3mm}
    \includegraphics[width=1.35in]{0411010hr0x4resized.png}
    \includegraphics[width=1.35in]{0411010hr1x4resized.png}
    \includegraphics[width=1.35in]{0411010hr2x4resized.png}
    \includegraphics[width=1.35in]{0411010hr3x4resized.png}
    \includegraphics[width=1.35in]{0411010hr4x4resized.png}\\
    \\
    \end{tabular}
\vspace{-5mm}

    \caption{\textbf{Comparison of \nj{output} images \nj{in} each stage.} \nj{For better comparison, each row shows a part \wacvadded{(64 $\times$ 64 pixels)} of an image} in DIV2K image set. \wacvadded{As described in the introduction, the originals are naturally blurred in the capturing and the transmission process.} The stage number increases \nj{from left to right}.  \nj{As the stage} increases, it can be observed that \nj{the texture is clearly seen and} the result becomes clearer. {It is easily identifiable in a man's beard, in the particle of the bricks, and in the clarity of the letters. }}
\label{result_in_image}
\vspace{-3mm}
\end{center}
\end{figure*}

\subsection{Setting the number of stages}
In our proposal, \nj{the number of stages $N$ can be set as the point where the quality of $HR_{N+1}$ is worse than that of $HR_{N}$ and the corresponding network $SR_N$ can be used as the final super resolution network.}
The important thing here is how to judge \nj{whether the output image} gets better or not. As a result of \nj{our} experiment, \nj{we can consider a couple of} solutions. The easiest approach is to measure the quality of the actual \nj{validation} set to see if this score is increased. This is intuitive, but it has the disadvantage that it is difficult to use it because of the time and effort required for measurement. The second method is to confirm whether the change \nj{$||HR_{N+1} - HR_{N}||$} becomes larger or smaller as the repetition progresses. As the output of each stage approaches a certain ideal point, the difference \nj{will} gradually decrease. When the output of the stage deviates from the ideal point, the divergence starts and the difference \nj{will increase}. This will be explained in more detail in the \nj{following experimental results}.

\section{Experimental result}

\subsection{Training setup}
\nj{In our experiment, we used the DIV2K image dataset \cite{DIV2K}, which contains 2K RGB images.
This is to prove \nj{the performance of image enhancement on contents in the FHD level which} occupies the majority of \nj{current} broadcasting \nj{contents}. The DIV-2K \nj{dataset} consists of 800 train sets and 100 validation sets. \nj{From the 800 DIV-2K training set, we made our \nj{training} database which was} constructed by ripping 100 random patches \nj{with a} 256$\times$256 size. After making \nj{the original images, $HR_0$}, we adjusted the JPEG compression and \nj{down-scaled the patches} to \nj{the} half size to create the LR images. \sam{The learning was carried out from the 80,000 pairs made in this way.}

\nj{
In all of our experiments, we used MSE as a loss function and Adam as an optimizer. Batch size of 8 was used with learning rate of 0.001. We trained the network for 20 epochs.
}

\subsection{Convergence of recurrent learning}
In order to confirm whether the \nj{proposed} recurrent learning \nj{converges} to the maximum improvement point or simply diverges as if the Laplacian filter is applied multiple times, we observed the \nj{difference} of \nj{output images between consecutive stages}. The \nj{image difference} is an objective indicator that shows the average difference of all the pixels in \nj{an} image in \nj{terms of ratio}. It can be similar to MSE, but it is normalized to the image's brightness and resolution. \nj{It is a more interpretable measure in that it is} more intuitive \nj{telling the percentage of the image}  changed. The average image difference of the \nj{100 validation images} in each stage \nj{is shown in} Table \ref{table:image_diff}.

\nj{Figure \ref{fig:image_diff_graph}} shows Table \ref{table:image_diff} to a graph. It shows that the learning converges to a certain point. And it can be seen that the \nj{change of difference starts to increase} after a specific point, the Stage 4.
In order to confirm that the recurrent learning is different from the repetitive application of \nj{a} simple filter, we \nj{repeatedly applied unsharp-masking (USM) filter to the original image, which is shown by} the green line. The difference \nj{increases} linearly and \nj{image difference} \nj{from Stage (N-1) does not converge ever}. Also, from the \nj{observation that the} image difference is linearly increased \nj{by repeated USM application}, we can confirm that the image difference is properly measured.

\begin{table*}
\caption{\textbf{Result of VIQET image quality measurement.} The final MOS score in this table is obtained based on the detailed score of 10 items excluding \textsc{Flat Region Index} which results in NaN for all the cases.
Except for over and under exposed ratio, a high score lead to a high MOS, which means a better quality. }
\begin{center}
\small
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
	&MOS&Resol&Edge&SN&Texture&Satu&Color&Illumi&Dynamic&Over&Under\\
	Name&Score&ution&Acutance&Index&Acutance&ration&Warmth&nation&Range&Expose&Expose\\
\hline\hline
    \texttt{LR}&	3.645&	2.8&	52.6&	311.9&	96.3&	116.8&	107.7&	156.8&	101.0&	0.1&	0.4\\
    \texttt{HR0}&	3.825&	2.8&	56.6&	314.8&	108.7&	117.3&	107.8&	163.6&	101.0&	0.1&	0.4\\
    \texttt{HR1}&	3.962&	2.8&	82.7&	322.3&	150.7&	117.7&	107.8&	182.7&	101.3&	0.1&	0.4\\
    \texttt{HR2}&	4.028&	2.8&	92.5&	325.3&	167.8&	\textbf{118.4}&	\textbf{108.0}&	190.0&	101.2&	0.0&	0.4\\
    \texttt{HR3}&\textbf{4.057}&2.8&\textbf{94.3}&	326.2&	\textbf{171.8}&	118.0&	107.8&	\textbf{191.4}&	101.2&	0.0&	0.4\\
    \texttt{HR4}&	4.050&	2.8&	93.1&	326.7&	171.4&	117.1&	107.2&	191.1&	101.3&	0.1&	0.4\\
    \texttt{HR5}&	4.048&	2.8&	92.1&	326.7&	170.5&	117.7&	107.6&	190.7&	101.3&	0.1&	0.4\\
    \texttt{HR6}&	4.038&	2.8&	89.5&	327.6&	168.1&	116.0&	107.8&	188.2&	101.1&	0.0&	0.4\\
    \texttt{HR7}&	4.030&	2.8&	87.9&	329.0&	166.9&	114.6&	107.3&	186.8&	101.1&	0.1&	0.4\\
    \texttt{HR8}&	4.030&	2.8&	87.3&	\textbf{331.1}&	166.1&	111.7&	106.1&	186.2&	101.2&	0.1&	0.4\\
\hline
\end{tabular}
\end{center}
\vspace{-5mm}
\label{table:mos}
\end{table*}

\subsection {The result images}
\nj{Figure \ref{result_in_image}} shows \nj{some samples of the output images in} each stage. It shows how the \nj{resultant images} for each stages are being improved. All images are down-scaled with the same resolution as the original, and we can visually confirm that the \nj{quality} of the images \nj{is} improved. The improvement level at each stage is similar to the slope of the graph in \nj{Fig. \ref{fig:image_diff_graph}}. \nj{As stage goes on}, it can be observed that the texture is clearly seen and the result becomes clearer.
{It can be said that recurrent learning converges to a certain level of improvement, since there is little difference between $HR_{3}$ and $HR_{4}$ \nj{in all the samples}.}

\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{mos.png}
\caption{\textbf{The MOS results for each stage.} The final score, MOS (Red line), peaked at $HR_{3}$. {The increase in MOS is similar to the convergence of the difference, yellow line in Fig.\ref{fig:image_diff_graph}.
Also the \nj{specific measures in other colors} tend to be increasing until $HR_{3}$ and then \nj{they saturate} or decrease slightly.}
}
\vspace{-3mm}
\label{fig:MOS}
\end{figure}

\subsection {Verification using VIQET tool}
In order to clearly demonstrate the effectiveness of image enhancement, we looked for ways to more objectively measure the image quality. Many non-reference quality measure models \nj{exists}, but most \nj{of them need} pre-learning the references and \nj{check} how close the measured image is from \nj{the} learned. This \nj{is} not what we want, because we \nj{do not} have any reference. Fortunately, the \textit{Video Quality Experts Group} (VQEG) has developed a tool that measures the image quality and \nj{developed} the Mean Opinion Score (MOS) of \nj{an} image as 0 (Poor) to 4.5 (Best). Using this tool, we can evaluate the image quality improvement for each stage objectively.

\nj{In the} VQEG image quality Evaluation Tool (VIQET) \cite{VIQET}, the average MOS score is derived for each image through a total of 11 measures:
{\small
\textsc{1. Resolution, 2. Multi-scale Edge Acutance, 3. Noise Signature Index, 4. Flat Region Index, 5. Multi-scale Texture Acutance, 6. Saturation, 7. Color Warmth, 8. Illumination, 9. Dynamic Range, 10. Over Exposed ratio, 11. Under Exposed ratio}}. \nj{To use VIQET, we must set a category of an image among 4 categories:
{\small
\textsc{Outdoor Day - Landscape, Indoor - Wall Hanging, Indoor - Arrangements}
} and
{\small \textsc{Outdoor Night - Landmark}}.
Most of the DIV-2K images \nj{are} bright images taken from the outside and we simply set all \nj{the images} to the \textsc{outdoor day} category.} The \nj{detailed average} scores of 100 DIV-2K \nj{validation} images for each stage are shown in Table \ref{table:mos}.

The key point is that the overall score, MOS value, increases with \sam{the stage}. The result proves that ``The image quality can be increased by recurrent learning''. The MOS results for each stage are shown in Fig. \ref{fig:MOS}. Based on the MOS, it can be confirmed that the image quality from the original \nj{3.825} was improved up to 4.057. The increase is 0.232 points \nj{from the original}, which is larger than the difference between LR and \nj{the} original, 0.180. \sam{It prove} that image enhancement through recurrent learning is efficient. \nj{There} is a limitation \nj{in improving the image quality based on the difference between $HR_0$ and LR with one-time learning and we have resolved this by providing a new target obtained through recurrently training the network.}
As can be seen in Fig. \ref{fig:MOS}, the MOS increase is largely due to the improvement of four terms. The first one is the improvement of multi-scale edge acutance \nj{value.} It means the edge got to be more sharp. The second is the Noise Signature Index, which means that the distinction between image and noise is more clear. The \nj{third} one is multi-scale texture acutance, and it means increment of {the level of activity and detail in the scene. Lastly, illumination means the light has become more enriched, also \nj{indicating} that the image has more detail.}

\begin{table}
\caption{\textbf{MOS results for various networks.} This table shows that there is an additional improvement in the index after HR1. It means the application of Recurrent Training Strategy (RTS) is meaningful regardless of the network. And we can confirm the proposed network is more suitable for recurrent learning. \nj{Images are available in the supplementary material.}}
\begin{center}
\small
\vspace{-3mm}
\begin{tabular}{|c||c|c|c|c|c|}
\hline
	Network&Param&HR0&HR1&HR2&HR3\\
\hline\hline
    \texttt{\textbf{Proposal}}&&&&&\\
	\texttt{Residual1}  & 41K   &2.384&2.639&2.701&2.712\\
	\texttt{Residual3}  & 114K  &2.384&2.637&2.719&2.720\\
	\texttt{Residual6}  & 225K&2.384&\textbf{2.644}&\textbf{2.725}&\textbf{2.744}\\
\hline
    \texttt{\textbf{Reference}}&&&&&\\
    \texttt{SRCNN}      &57K    &2.384&2.603&2.680&2.680\\
    \texttt{VDSR}       &666K   &2.384&2.599&2.668&2.669\\
    \texttt{SRResNet}   &1.5M   &2.384&2.617&2.720&2.703\\
    \texttt{EDSR(64$\times$16)}&1.7M   & 2.384 &2.632&2.725&2.737\\
\hline
\end{tabular}
\end{center}
\vspace{-7mm}
\label{table:comparison}
\end{table}

\iffalse
\begin{figure}
\includegraphics[width=\linewidth]{comparison_test.png}
\caption{\textbf{The MOS results for various networks.} Although the rising speed of MOS is different for each algorithm, the tendency of learning is similar. The difference between HR1 and HR2 is bigger than the difference between networks. {It means that it is more dependent on the number of learning times than the size of the network, and as a result, it can produce good results even with light networks.}}
\label{fig:compasion_network}
\end{figure}
\fi

\subsection{RTS \nj{on} various networks}

\nj{To investigate the characteristics of the proposed RTS more clearly and to check whether the proposed RTS} is still effective with other \nj{algorithms}, we test \nj{RTS on} several other networks
as well as the proposed network.
\nj{RTS was applied to total} six networks \nj{and the performances are shown in Table \ref{table:comparison}}. \nj{In the table, the above three networks are designed based on our proposal shown in Fig. \ref{fig:network} by changing the number of residual layers in the middle to 1, 3 and 6, respectively.}
The proposed architecture in Fig. \ref{fig:network} has the advantage of being able to easily compare the number of parameters in the network \nj{considering the number of repetitions}.  For example, a comparison between learning once with doubled parameters and learning twice with halved parameters \nj{is possible}. In the three bottom rows,  we have \nj{chosen} SRCNN, VSDR, \wacvadded{SRResNet and ESDR} to check whether our \nj{recurrent training strategy} is meaningful even in \nj{a} heavier super-resolution network. \nj{By this experiment, the characteristics} of \nj{the proposed RTS on various networks can be checked} from the simplest network, SRCNN, to \wacvadded{a} latest heavy network, \wacvadded{EDSR}.

Unfortunately, the existing networks are so heavy that the GPU memory \nj{overflows}, making it impossible to keep the original size of a DIV2K image intact. So we have to divide each image into 4$\times$2 and processed \nj{them independently}. The MOS result is \nj{the} averaged value of 800 pieces of divided 100 \nj{validation images}. The segmentation of the image caused a decrease in the content and resolution, which resulted in a decrease in the MOS score. However, since this segmentation is applied equally to all the six \nj{networks}, the tendency is comparable and the experiment is fair. Table \ref{table:comparison} shows the resultant MOS scores as well as the number of parameters in each network.

In \nj{Table} \ref{table:comparison}, the scores of HR2 and HR3 are higher than \nj{HR1 at all times}. This suggests that recurrent learning may lead to additional performance improvement in existing algorithms. In addition, it is shown that the improvement \nj{due from} the learning method is more meaningful than the improvement \nj{caused by using a different} network. It means that the \nj{quality of an image} is more dependent on the number of repetitions in RTS than the size of the network, and as a result, the proposed RTS can produce good results even with light networks.
Although the \nj{amount of increase in} MOS is different for each algorithm, the tendency of learning is similar. From this, we can say that applying \nj{RTS} makes it possible to use \nj{a} less expensive network.

\section{Conclusion}

In this paper, we proposed a new image enhancement method by a recurrently-trained super-resolution network.
\nj{The proposed recurrent learning} consists of two phases: \nj{Network Training} phase to learn HR form LR, and \nj{Target Update} phase to apply \nj{the trained network} to the original to produce HR \nj{image for} the next stage. We clarified the \nj{characteristics} of recurrent learning by analyzing the process and results. Using image difference, we can show the results converge to a specific improvement point without divergence. Also, by using VIQET mean-opinion-score, it has been \nj{numerically shown that the quality of an image improves clearly as we repeat the learning procedure.}
It is the first time to use these objective \nj{measures} in the image-enhancement \nj{area}, \nj{which is meaningful in} that the MOS score \nj{is useful} not only in a theory \cite{1248017} but also in practice.

\nj{The DIV2K image dataset used in this paper has an equivalent image quality to that of an} actual digital broadcasting video, \nj{which is different from the conventional datasets with low quality used in the existing super-resolution papers.} Nonetheless, by RTS, the MOS score increased more \nj{than} the difference between the original and its low-resolution version, which is impossible with the conventional one-time learning. The increased MOS score proves that image enhancement through recurrent learning is efficient.

In order to clarify the \nj{effectiveness of the proposed RTS}, we experimented with various networks. The results show that recurrent learning can make an additional improvement, regardless of the  network used. And the results show a way to \nj{obtain a} cheap but effective network. \nj{We expect that it will contribute to the related society by dramatically reducing the cost of super-resolution and image enhancement.}

{\small
\bibliographystyle{ieee}


\section{Introduction}
\label{sec:intro}
Computer simulators \cite{schonlau1997computer} and/or sophisticated in-house laboratory experiments \cite{flournoy1993clinical} representing physical phenomena generally pose the problem of being computationally or logistically challenging. In scenarios where extracting information \cite{jaynes1957information,kullback1997information} about an underlying physical process to forecast or predict its behaviour is critical, generating hundreds of thousands of experimental data points becomes infeasible. To alleviate the above issue, surrogate modeling (also called predictive modeling) is a popular approach towards obtaining an inexpensive representation of the simulation or experiment.

Among the myriad surrogate modeling \cite{gunst1996response} techniques, Gaussian Process (GP) \cite{williams2006gaussian} regression (GPR) is a well established technique in the area of probabilistic modeling of expensive functions. GPR is a non-parametric machine learning technique in the sense that it does not make assumptions about the functional form of the function being modeled. However, the GP model does have parameters, called hyperparameters, that capture characteristics such as periodicity, smoothness, measurement noise, etc., of the underlying function. These hyperparameters are learned from the data (limited set of observations) using empirical Bayesian techniques like maximum likelihood estimation
\cite{pan2002maximum}, or fully Bayesian \cite{kennedy2001bayesian} methods based on Markov chain Monte Carlo (MCMC) \cite{robert2013monte,gelman1995bayesian}. For a typical GP the number of hyperparameters is on the order of 10-100.
Furthermore, the predictive accuracy of the GP is highly sensitive to the values of the hyperparameters thus requiring sophisticated optimization. The high dimensionality of the hyperparameter space prompts a rigorous MCMC treatment \cite{cowles1996markov} to learn the hyperparameters. Also, the Bayesian MCMC \cite{gilks1995markov,carlo2004markov,green1995reversible} methods allow quantification of uncertainties \cite{ghahramani2003bayesian} around the estimates of the hyperparameters, which makes possible drawing samples from the space of the underlying function.

Recent work on extending GPs to big-data applications has focused on deriving variational representations of GPs~\cite{peng2017asynchronous}, constructing sparse approximations of GPs~\cite{hensman2013gaussian,quinonero2005unifying} and training local GPs using informative subsets of the data~\cite{snelson2007local,lee2017hierarchically}.
However, for most of these methods, the modeling assumes the parameters of the covariance kernel fixed to point estimates that are obtain using optimization or are considered known.
A recent review of methods for extending GPs for big-data applications can be found in~\cite{liu2018gaussian}.
In this paper, we restrict our focus to Bayesian GPs where the covariance parameters are inferred using MCMC.

GE's Global Research Center has its own implementation of fully Bayesian modeling, called GE Bayesian Hybrid Modeling (GEBHM), which has been applied \cite{ghosh2018bayesian,kristensen2016expected} successfully to various engineering problems over many years.
Given the greater predictability achieved by using MCMC, the  MCMC method comes with its own pitfalls \cite{o1987monte}.
For example, in problems involving huge training data sets, on the order of $N=1000$ data points, the excessive computational training time is a direct consequence of a large number of GP covariance matrix inversion operations, scaling as $O(N^{3})$. The work done in \cite{leithead2007n} has focused on circumventing this problem by obtaining approximations to the inverse of the matrix. Scalability of the MCMC methods with varying complexity of the problem, such as increasing number of input dimensions \cite{haario2004markov}, multiple correlated functions, etc., can break down when applied to challenging industrial problems.

 Sequential Monte Carlo (SMC) \cite{doucet2001introduction,diaconis2003sequential} methods offer a promising alternative considering the availability of multiple individual processor elements (PEs), also called cores, and high performance computing environments. The inherent parallelizability paves the way for SMC methods to be used in fitting the hyperparameters through a form of importance sampling and reduce the required number of matrix inversions per PE. Moreover, SMC methods can be tweaked by the user based on the type of problem through multiple \emph{knobs} which determine characteristics of the SMC such as number of simultaneously running MCMC chains, number of MCMC steps per chain, discretization of the posterior distribution into finite samples, etc. This magnanimity is an innate feature of most SMC algorithms \cite{andrieu2001sequential,andrieu2010particle}, the theory of which guarantees convergence under specific assumptions \cite{bilionis2012free,bilionis2015crop}. We apply an adaptive SMC (ASMC) method \cite{kristensen2013relative}, by leveraging \emph{pySMC} (a package written in the PYTHON programming language). The details of the pySMC library are provided at \url{https://github.com/
ebilionis/pysmc}. The ASMC method samples according to the same kernel as GEBHM's standard MCMC, albeit reducing the number of covariance matrix inverse computations by an order of magnitude. This difference is critical to the ASMC alleviating the burden of excessive computational time in GP model training. We extend the capabilities and flexibility of the ASMC by allowing for a predefined set of importance sampling distributions which could be provided by the user, based on the complexity of the problem. As an addendum to the methodology, fine-tuning of the widths of proposal distributions of the hyperparameters is done on the fly, independently.

The purpose of this study is to verify the scalability of ASMC, its higher predictive accuracy, and to demonstrate the savings in time (which translate to cost savings) that can be leveraged henceforth. The investigation is carried out using six problems with varying types of complexity in terms of the input dimensionality, the number of outputs, the size of the training data, etc. The comparison tests are done using multiple PEs on a workstation computer to highlight the ease of implementation and usability of the ASMC. Further comparison tests include the results obtained by the seamless scalability achieved by ASMC using high performance computing (HPC) clusters. Out of the six chosen problems four are synthetic in nature with different dimensionality and training data size and two stem from GE industrial applications.

The outline of the paper is as follows. We start \sref{metho} by providing some basic mathematical definitions of GPR and the details of the posterior of the
hyperparameters of the GP. Our numerical results are presented in \sref{results}. The impact of the ASMC on problems with large training data sets is highlighted in \sref{scalability}.
In particular, in \sref{toy_small}, \sref{toy_medium} and \sref{toy_large}, we validate our approach using three synthetic
problems with varying input dimensionality.
In \sref{flame_transfer} and \sref{combustions}, we apply the ASMC methodology to solve a challenging steam turbine compressor problem and a multi-objective combustion
problem, respectively. We present our conclusions in \sref{conc}.

\section{Methodology}
\label{sec:metho}
\subsection{Gaussian process regression}
\label{sec:gp}
The surrogate models that GEBHM builds for the overall objective(s) for the problem are data-driven statistical models that rely on the Bayesian nature of GPR \cite{williams2006gaussian}.
GPs are non-parametric meaning, they do not assume a model form (think: polynomial of degree $n$, among others). GPs try to map inputs to outputs by superposing several Gaussian distribution functions. GPs can accurately model complex variations in the data, they can capture discontinuities effectively and can produce very accurate surrogate statistical models with comparatively less data. An important feature is that GPs offer a measure of the quality of the surrogate statistical model and direct us to areas of input space where we need more data, for example \cite{kristensen2017adaptive}. We now go over the finer details of GPR where we define the posterior distribution of the hyperparameters which we want to sample from.
We use a zero mean for the GP. The covariance of the GP is computed based on a kernel that is a function of the distance between the inputs. For a multi-objective problem the covariance matrix is a block diagonal matrix comprising of the individual covariance matrices for each of the objectives. The covariance for an $m$ output, $d$ input, problem with $N$ training data points is shown below:
\begin{equation}
\Sigma=
\begin{bmatrix}
\Sigma^{(1)} &		&		\\
				& \ddots & 	\\
				&		&\Sigma^{(m)}
\end{bmatrix},
\label{eqn:cov_overall}
\end{equation}
where each of the individual elements along the diagonal represent the covariance matrix of dimension $N \times N$ of a single output.
These can be more generally expressed as in \qref{cov_ind}.
\begin{equation}
\Sigma^{(k)}_{ij} = \frac{1}{\lambda_{z}^{(k)}}\sum_{l=1}^{d} \exp\left[-\beta_{l}^{(k)}{(x_{i,l}-x_{j,l})^2}\right] + I\frac{1}{\lambda_{s}^{(k)}} + I\frac{1}{\lambda_{o}}
\label{eqn:cov_ind}
\end{equation}
where $i,j \in \{1,\cdots, N\}$, $ k \in \{1, \cdots, m\}$, $l \in \{1, \cdots, d\}$, $x_{i,l}$ is the value of the $i$th training point in the $l$th input dimension, and let $\bm{\beta}^{(k)}=\{\beta_{1}^{(k)} \cdots \beta_{d}^{(k)}\}$.
The hyperparameters of the GP model are $\bm{\theta} = \{\bm{\beta}^{(k)}, \lambda_{z}^{(k)}, \lambda_{s}^{(k)},  \lambda_{o}\}_{k=1}^{m}$.
The $\bm{\beta}^{(k)}$ parameters represent the inverted length-scales of the covariance kernel. $\lambda_{z}^{(k)} $ is also called the signal-strength which is proportional to the scale of the magnitude for the objective space and $\lambda_{s}^{(k)}$ represents the precision of the random Gaussian measurement noise. $\lambda_{o}$ represents the precision of the random Gaussian measurement noise to account for the overall randomness between the multiple outputs. It can be inferred from above that for our problems of learning the model's hyperparameters, it is always a problem of estimating $m\times (d+2) + 1$ number of parameters.
The likelihood of observing the training data for a selected set $\bm{\theta}$ of hyperparameters is:
\begin{equation}
L(\bm{Y}|\bm{\theta}) = \frac{1}{|\Sigma|^{\frac{1}{2}}}\exp(-\frac{1}{2}\bm{Y}^{T}\Sigma^{-1}\bm{Y}),
\label{eqn:likelihood}
\end{equation}
where the $i, j$th element of the matrix $\bm{Y}$ is the $j$th output value for the $i$th training data point.
The conditional posterior of the hyperparameter set $\bm{\theta}$ can be written as follows:
\begin{equation}
p(\bm{\theta}|\bm{Y}) \propto L(\bm{Y}|\bm{\theta})p(\lambda_{o})\prod_{k=1}^{m}p(\bm{\beta}^{(k)})p(\lambda_{z}^{(k)})p(\lambda_{s}^{(k)}).
\label{eqn:posterior}
\end{equation}
The posterior given in \qref{posterior} is the so-called \emph{target distribution} known only up to a proportionality constant. The sequence of sampling distributions in the ASMC are derived from the functional form given in \qref{posterior}.
\subsection{MCMC in GEBHM}
\label{sec:bhm_mcmc}
The implementation of MCMC in GEBHM works in two phases: a) the initialization phase and b) the main chain.
The initialization phase of the MCMC solves the purpose of selecting the proposal width of the proposal distributions for the main (full) chain of the MCMC. This initialization process can only be performed on a single PE, thus limiting the computational power offered by multiple cores. The main chain can use multiple cores, but the single core operation of initialization does contribute to a significant percentage of the total cost of the MCMC in GEBHM. The complete details of the initialization process have been omitted to protect company proprietary information. The MCMC algorithm uses a Metropolis algorithm \cite{gelman1996efficient} for making jumps (steps) during the main chain.
\subsection{Adaptive sequential Monte Carlo}
\label{sec:asmc_math}
The sequential Monte Carlo methodology (also known as particle filtering in some literature) \cite{gordon1993novel} essentially uses a coherent sequence of distributions to do MCMC sampling from. This sequence of distributions starts with a distribution mimicking a uniform distribution and progressively moves closer (in analytical form) to the posterior distribution \qref{posterior} (in our case).
A sequence of $n$ sampling distributions based on the empirical form given in \qref{pi_gamma} can be built with $\gamma_{0} \sim 0 $ and $\gamma_{n} = 1$. Thus, the sampling process can be understood as a slow progression from a non-informative \emph{uniform} distribution to the target posterior distribution  \qref{posterior}. The jumps between states for all sampling distributions are based on a Metropolis stepper (same as the one being used for MCMC in GEBHM).
The conditional posterior in \qref{posterior} is reproduced in \qref{smc_true} while we shun the dependence on $\bm{Y}$ to clearly represent the form of the sampling distributions of the ASMC as a function of only the hyperparameters $\bm{\theta}$.
The sampling distributions are illustrated in \qref{smc_dist}. Each $p_{i}(\bm{\theta})$  is a sampling distribution and each chain (ASMC particle) takes $num_{asmc}$ steps before moving on to a different value of $\gamma$.
\begin{eqnarray}
p(\bm{\theta})&\propto&\pi(\bm{\theta})
\label{eqn:smc_true}\\
{p_{i}(\bm{\theta})}&\propto&{\pi_{i}(\bm{\theta})}
\label{eqn:smc_dist}\\
{\pi_{i} (\bm{\theta})}&=&{\pi(\bm{\theta})}^{\gamma_{i}}.
\label{eqn:pi_gamma}
\end{eqnarray}
At each $\gamma_{i}$ the corresponding distribution ${p_{i}(\bm{\theta})}$ is described as a discrete particle approximation as shown in \qref{particle_approx}.
\begin{equation}
{p_{i}(\bm{\theta})} \approx {\mathlarger{\sum}}_{j=1}^{N}w^{j_{i}}\delta(\bm{\theta}-\bm{\theta}^{j_{i}})
\label{eqn:particle_approx}
\end{equation}
These finite number of samples are called particles (parallel MCMC chains). At each value of $\gamma$ each particle takes a predefined number of steps $num_{asmc}$. Thus, each particle evolves its own MCMC chain based on the Metropolis stepper.
The values of $\gamma$ visited by the algorithm are governed by a statistic known as the \emph{effective sample size} (\emph{ess}). The \emph{ess} is calculated as defined in \qref{ess}. Thus, when all particles have equal weights the \emph{ess} is equal to the number of particles, $N$.  The \emph{ess} guides the methodology finding the next $\gamma$ i.e. $\gamma_{i+1}$  as shown in \qref{ess}.
\begin{eqnarray}
ess|\gamma_{i}&=&{\mathlarger{\sum}}_{j=1}^{N}\frac{1}{{w^{j_{t}}}^{2}}
\label{eqn:ess}\\
ess|\gamma_{i+1}&=&ess\_reduction \times ess|\gamma_{t},
\label{eqn:ess_reduction}
\end{eqnarray}
where $ess_{reduction}$ is a predefined threshold $\in{(0, 1)}$.
The $\gamma_{i+1}$ is usually selected via an algorithm like \cite{brent1980improved}.
For more details on the \emph{ess} and the selection of $\gamma_{i+1}$ the reader is referred to the work of \cite{bilionis2015crop,bilionis2012free}.
The dynamic selection of the set of $\gamma$s for a single ASMC requires special attention in the industrial context.
For some problems the dynamic selection process can be too time consuming.
Towards ameliorating this we replace the dynamic selection of the $\gamma$ parameter with a prespecified set of $\gamma$s called the \emph{grid}. The \emph{grid} used in this works allows uniform sized jumps for the $\gamma$ parameter but we expect better results with a non-uniform \emph{grid}.
Consider \qref{ess_reduction}, the number of $\gamma$s visited by the ASMC can vary based on factors such as the size of training data, number of input dimensions, etc. The \emph{grid} thus constrains the freedom of the algorithm especially when time is a greater priority for a user than small gains in accuracy.
We list the steps taken by GEBHM's ASMC methodology in algorithm \ref{alg:asmc}.
\begin{algorithm}[htb]
\caption{ASMC in GEBHM.}
\begin{algorithmic}[1]
\Require
		   number of steps per particle $num_{asmc}$,
		   \emph{grid}: $\bm{\gamma}$,
		   MCMC stepper,
	 \State Initialize the particles (chains) via MCMC stepper.
            \For {$i \in [1, n]$}
                \State Do $num_{asmc}$ steps for each particle.
                \State \emph{tune} the proposal widths of the proposal distributions for each hyperparameter based on the			respective acceptance ratios from the MCMC stepper.
                \State $\gamma \leftarrow \bm{\gamma}_{i}$
    			\EndFor
\end{algorithmic}
\label{alg:asmc}
\end{algorithm}
The added advantage of the ASMC is the lack of need to pre-compute the proposal widths (as is the case for MCMC in GEBHM) as it adaptively scales the proposal widths on the fly for each parameter. The exploration of posteriors of parameters with multiple modes is a problem that the ASMC is semantically designed to tackle.
\section{Results}
\label{sec:results}
We first present our numerical results demonstrating the ASMC's ability to scale well with large training data sets.
\subsection{Toy problem 1: Scalability}
\label{sec:scalability}
As mentioned in \sref{intro}, the standard MCMC implementation makes the GP model training process cumbersome. In cases with large training data, the aforementioned problem worsens with the whole data set not being put to use for training the GP model. Reducing the frequency of the computationally challenging part is critical to the consequence of applying a MCMC method to GP model training. We apply the ASMC to the following 10 dimensional non-linear mathematical problem:
\begin{eqnarray}
y(\mathbf{x})&=&3\sin(x_{1})x_{2} + \cos(x_{3})\sin(x_{4}) + \sin(x_{5})\sin(x_{6}) \nonumber \\
&& + \sin(x_{7}) + \sin(x_{8}) + 7x_{9} + 6x_{10}
\label{eqn:scalability}
\end{eqnarray}
\qref{scalability} enables us to generate a large training data set (greater than 1000 points).
The performance of the two algorithms, GEBHM and ASMC, is highlighted in \fref{scalability} (a),  (b) and (c). ASMC's ability to
handle large sets of data is clearly visible in \fref{scalability} (a) where the standard MCMC takes almost four times as long as the ASMC with around 1000 training data points. Secondly, the ASMC's GP model building times are comparatively low even with training data size greater than 3000 points. This is accomplished due to ASMC's seamless operational scalability to HPC environments. \fref{scalability} (b) highlights the predictive accuracy achieved by ASMC (which is mainly a satisficing objective) being equal to or better than the standard MCMC.
The MCMC settings were kept at 6800 steps including 2000 steps for the initialization phase.
The savings in absolute time shown in \fref{scalability}(c) provide an idea of the potential improvement made by using ASMC for GP model building with high training data problems.
It is important to consider the non-parallelizability of the initialization (shown in purple) of the MCMC, which limits the depth of impact that can be achieved by using a larger number of cores (greater than 12) for the GEBHM. The times represented by the initialization will remain fixed irrespective of the number of cores, thus making the comparison with ASMC even more realistic.
\begin{figure}[htp]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{time_data_scalability}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_data_scalability}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{time_savings_data_scalability}.png}}}
\caption{Subfigures~ (a), (b) and (c) show the scalability, predictive accuracy and savings in time for the ASMC compared to MCMC for a large training data set problem respectively.  The number of particles for the ASMC is  480, for all runs.}
\label{fig:scalability}
\end{figure}
For the scalability problem we allow the ASMC to identify the set of sampling distributions on the fly, see \qref{ess_reduction}, hence not specifying the \emph{grid}.
\subsection{ASMC and MCMC settings}
Now, we move on to the comparison of the ASMC with GEBHM's MCMC on three synthetic problems. We plot the RMSE for the ASMC obtained using the final particle approximation for different number of particles, keeping the number of cores on a workstation fixed. As the number of particles increases, the time taken to build the GP model also increases and the corresponding RMSE shows gradual decrease. This is done because plotting the RMSE by using particles from the different sampling distributions of a single ASMC run would be misleading simply because the sampling distributions visited by the ASMC are different from one another. This is a crucial step, where the ASMC differs from the conventional MCMC (where the sampling distribution remains the same i.e. the posterior distribution of the hyperparameters known up to a proportionality constant \qref{posterior}). \\
For workstation computations the ASMC settings are fixed to 6 cores, and the number of particles at 6, 12, 30, and 60. For HPC computations we vary the number of cores and fix the number of particles equal to the number of cores for each run of the ASMC. The number of cores for the HPC runs of the ASMC for all problems are 24, 48, 120, 240, and 480. For ASMCs performance on the HPC with greater than 480 particles (for the toy problem in \sref{toy_large}), the number of cores is fixed at 480 (the highest number available during our experiments on the HPC).

However, for the standard MCMC, the RMSE on the test data is plotted against time as the chains evolve. The settings for the standard MCMC are 5800 steps including 1000 steps for initialization. The $\gamma_{0}$ values are 0.001 and 0.00001 for the ASMC's runs on the workstation and the HPC, respectively.
\subsection{Toy problem 2}
\label{sec:toy_small}
The four dimensional problem, is a single-objective problem with a known analytic form of the objective. The objective function is a quadratic function of the four input variables. For safeguarding proprietary information of GE's business partners, the full functional form of the input-output relation has been omitted from the text.\\
\begin{figure}[h]
\centering
{\includegraphics[width=0.5\textwidth]{{{rmse_time_small_pc}.png}}}
\caption{Root mean squared error versus time taken to build the GP model for the two algorithms for the four dimensional toy problem. Number of particles (red dots) are 6, 12, 30, and 60, across the time axis for ASMC.}
\label{fig:toy_small}
\end{figure}
The purpose of comparing the two algorithms on a low-dimensional problem is to demonstrate the ability of the ASMC to replace the standard MCMC on a workstation, before showing comparisons between the two methodologies by leveraging HPC platforms.
\fref{toy_small} shows the savings in time on a workstation with the ASMC using 6 cores and the GEBHM standard MCMC using 12 cores after the completion of the initialization (done on a single core). Each red dot in \fref{toy_small} represents the RMSE on test data after the completion of the ASMC GP model building, as the number of particles for the ASMC is varied across time.
For the standard MCMC, the RMSE on test data during initialization is separated from the rest of the particles by the vertical dotted line. The standard MCMC RMSE versus time is plotted by computing the RMSE at each time step of the MCMC chain.
Higher dimensional problems are used to compare the two algorithms with ASMC scaled-up on the HPC.
The number of steps per ASMC particle is equal to one. The \emph{grid} is fixed at 10 steps and 20 steps for the workstation and HPC, respectively.
\subsection{Toy problem 3}
\label{sec:toy_medium}
We use an 18 dimensional, single-objective problem to further demonstrate the effectiveness of the methodology. The problem is a torsion vibration problem where the high natural frequency of a three-shaft and two-disk system is the objective function. The input variables include the diameter, length, rigidity and weight density of the three shafts and the diameter, thickness, weight density of the two disks. A diagram representing the structure of the system is shown in \fref{toy_medium_pic}.
The problem can be mathematically represented as follows:
\begin{eqnarray}
y&=&\frac{\sqrt{\frac{-b + \sqrt{b^{2} -4ac}}{2a}}}{2\pi}  \nonumber \\
a&=&1 \nonumber \\
b&=&-\Bigg(\frac{K_{1} + K_{2}}{J_{1}} + \frac{K_{2} + K_{3}}{J_{2}}\Bigg) \nonumber \\
c&=&\frac{K_{1}K_{2} + K_{2}K_{3}+K_{3}K_{1} }{J_{1}J_{2}} \nonumber
\label{eqn:toy_medium}
\end{eqnarray}
where,
\begin{eqnarray}
K_{i}&=&\frac{\pi G_{i}d_{i}}{32L_{i}} \nonumber \\
M_{j}&=&\frac{\pi t_{j}\rho_{j}D_{j}}{4g}  \nonumber \\
J_{j}&=&\frac{1}{2}M_{j}\Bigg(\frac{D_{j}}{2}\Bigg)^{2}
\label{eqn:toy_medium_params}
\end{eqnarray}
for $i \in \{1, 2, 3\}$ and $j \in \{1, 2\}$. The $d_{i}s$, $L_{i}s$, $G_{i}s$ and $\lambda_{i}s$ represent the diameter, length, rigidity and the weight density of the three shafts and $D_{j}s$, $t_{j}s$ and $\rho_{j}s$ represent the diameter, thickness, and weight density of the two disks, respectively.\\
The convergence analysis of the two algorithms is shown in \fref{toy_medium}.
The number of steps per ASMC particle is five. The \emph{grid} is fixed at 10 steps and 20 steps for the workstation and HPC, respectively.
\begin{figure}[h]
\centering
{\includegraphics[width=0.5\textwidth]{{{toy_medium_pic}.png}}}
\caption{The setup of the torsion vibration problem.}
\label{fig:toy_medium_pic}
\end{figure}
\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_medium_pc}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_medium_hpc}.png}}}
\caption{Subfigure~(a) Number of particles on a workstation (red dots) are 6, 12, 30, and 60, and Subfigure~(b) number of particles on the HPC (red dots) are 24, 48, 120, 240, and 480, across the time axis, respectively for the torsion problem.}
\label{fig:toy_medium}
\end{figure}
\subsection{Toy problem 4}
\label{sec:toy_large}
A 100 dimensional non-linear single-objective problem poses a strong challenge in terms of the number of input dimensions to test the performance of both methodologies. The explicit details of the functional relationship have been omitted from the manuscript.\\
The convergence is terms of the RMSE on test data as a function of time is shown in \fref{toy_large}. The convergence for the workstation based runs can be seen in \fref{toy_large} (a), where the ASMC reaches almost similar accuracy as the standard MCMC in less than half the time (with 30 particles). Similar convergence trends can be seen for the ASMCs runs on the HPC with varying number of particles (24, 48, 240, 480, 960 and 2400) with time.
The number of steps taken by each ASMC particle is one for both (workstation and HPC) cases. The \emph{grid} is fixed at 10 steps and 20 steps for the workstation and HPC, respectively. The finer \emph{grid} considered for the HPC is to gain greater predictive accuracy while relying on the HPC (large number of cores) to compensate for the computational overheads incurred.
\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_large_pc}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_large_hpc}.png}}}
\caption{Subfigure~(a) Number of particles on a workstation (red dots) are 6, 12, 30, and 60, and Subfigure~(b) number of particles on the HPC (red dots) are 48, 240, 480, 960 and 2400, across the time axis, respectively for the 100 dimensional problem.}
\label{fig:toy_large}
\end{figure}
\section{Industry problems}
\subsection{Steam turbine compressor problem}
\label{sec:flame_transfer}
This is a five dimensional, single-objective design optimization problem, where the training data set is observed from temperature experiments. The goal is to obtain the best system design given requirements on the temperature distributions.
The number of steps taken by the ASMC particles during the workstation based runs is five and for the runs on the HPC is one. The \emph{grid} is fixed at 10 steps for both, workstation and HPC.
\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_flame_transfer_pc}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_flame_transfer_hpc}.png}}}
\caption{Subfigure~(a) Number of particles on a workstation (red dots) are 6, 12, 30, and 60, and Subfigure~(b) number of particles on the HPC (red dots) are 24, 48, 120, 240, and 480, across the time axis, respectively for the steam turbine compressor problem.}
\label{fig:flame_transfer}
\end{figure}
The savings in time and, for this case noticeable, gains in predictive accuracy are shown in \fref{flame_transfer}. Not only does the ASMC show explicit time savings, but with the HPC  it also results in greater predictive accuracy than the GEBHM standard MCMC as shown in \fref{flame_transfer} (b). This indicates the dense representation of the sampling distributions made possible by using the HPC (where each particle runs on a separate core).
\subsection{Combustion problem}
\label{sec:combustions}
The problem is a GE combustion test data problem, where we aim to build a GP model to predict two emission quantities as a function of three measured temperatures, three fuel flow parameters, and air flow. Thus, there are 9 hyperparameters that need to be estimated for each individual objective. The total number of hyperparameters for the GP model is 19. The MCMC settings for this problem are: 10000 steps with 1000 steps during the initialization.
\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_combustions_y1_pc}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_combustions_y1_hpc}.png}}}
\caption{Subfigure~(a) Number of particles on a workstation (red dots) are 6, 12, 30, and 60, and Subfigure~(b) number of particles on the HPC (red dots) are 24, 48, 120, 240, and 480, across the time axis, respectively for the first objective of the combustion problem.}
\label{fig:combustions_1}
\end{figure}
\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_combustions_y2_pc}.png}}}
\subfigure[]{\includegraphics[width=0.5\textwidth]{{{rmse_time_combustions_y2_hpc}.png}}}
\caption{Subfigure~(a) Number of particles on a workstation (red dots) are 6, 12, 30, and 60, and Subfigure~(b) number of particles on the HPC (red dots) are 24, 48, 120, 240, and 480, across the time axis, respectively for the second objective of the combustion problem.}
\label{fig:combustions_2}
\end{figure}
The significance of the results in \fref{combustions_1} and \fref{combustions_2} is noticeable from the perspective of the number of objectives. Not only does the scaled-up HPC ASMC do well in terms of computational time but it improves the predictive accuracy of the GP model compared to GEBHM's MCMC. The number of steps per ASMC particle is one for the workstation and five for the HPC runs. The \emph{grid} is fixed at 10 steps for both computing environments.

This provides a holistic overview of ASMC's performance on different types of problems. Problems of varying input dimensionality, varying number of objectives, and different training data sizes have been treated with the extended ASMC and compared to the current MCMC implementation. The ASMC does equally well or better compared to the MCMC in terms of predictive accuracy, while saving time by half or more in most of the challenging problems.
\section{Conclusions}
\label{sec:conc}
We demonstrate the working of the implementation of an ASMC method in GEBHM, and compare it to the performance of the existing semi-parallelizable MCMC algorithm on mathematical problems and challenging industrial problems with varying complexity. The main reason behind leveraging ASMC is the inherent flexibility and scalability offered by the algorithm's parallelizable nature, which is practically implementable with the availability of \emph{high-end} computational resources such as HPCs. The numerical results on toy problems of different input dimensionality show the versatility of the extended ASMC on both, workstation and HPC environments. This is augmented by extending the flexibility of the ASMC to allow the user to prespecify the sequence of sampling distributions. Thus, allowing the algorithm to balance savings in time and predictive accuracy. It is interesting to note that the ASMC has the inherent capability to treat different problems differently, but can also almost always perform well with a reasonable setting of the default parameters. The different parameters chosen with different problems in \sref{results} further highlights this nuance of the ASMC.
ASMC's performance on the combustion data problem, demonstrates its ability to perform equally well with multiple objectives.
Future directions of work on studying the algorithm include determining an optimal set of parameter settings for an arbitrary problem. This is important from the perspective of making the ASMC's implementation deliverable to users averse to the flexibility of the ASMC.
\section{Acknowledgements}
The authors wholeheartedly thank Dr. Ilias Bilionis, Purdue University for providing technical insights into this research.
\nolinenumbers