\begin{figure*}[htb!]
\includegraphics[width=1.0\textwidth]{figures/workflow.png}
\caption{In the first step, the user is presented with a gallery and scores the shown materials according to their taste. Then, a regression is performed to obtain a preference function via Gaussian Process Regression, which can be efficiently sampled for arbitrarily many new material recommendations. These recommendations can be visualized in real time using our neural network in a way that closely resembles the images rendered with global illumination. In the final step, the recommended materials can be conveniently assigned to an existing scene.}
\label{fig:workflow}
\end{figure*}

\section{Overview}
\label{sec:realoverview}

The overall workflow of our system consists of two stages (see also Fig.~\ref{fig:workflow}). In the \emph{first stage}, the user is presented with a gallery and asked to assign scores to the shown materials. After choosing a threshold value to control the variability of the output, a set of recommendations are computed (Section \ref{sec:gpr}) and visualized via neural rendering (Section \ref{sec:cnn}). \Revision{The recommendations depend entirely on the user scores and may span multiple material classes. If the user wishes to fine-tune a subset of the recommended materials,} a latent space is inferred for low-dimensional exploration (Section \ref{sec:gplvm}). The user then evaluates the result:
%\vspace{-5mm}
\begin{enumerate}
	\item If the recommendations are acceptable -- proceed to the next stage,
	\item If the recommendations need refinement -- assign scores to the \emph{newly proposed gallery} or \emph{adjust past rankings} and compute a new round of recommendations.
%	\item If the recommendations are not acceptable - re-assign the \emph{scores} and compute a new round of recommendations.
\end{enumerate}
%\vspace{-5mm}
In the \emph{second stage}, the user can choose from two ways to assign the recommended materials to a scene:
%\vspace{-5mm}
\begin{enumerate}
	\item Automatic workflow -- randomly assign the materials to the selected objects in the scene (Figures \ref{fig:teaser} (right) and \ref{fig:wrecking}). This is ideal for mass-scale material synthesis, when hundreds of materials are sought,
	\item Assisted workflow -- assign the materials to the scene manually and perform fine-tuning via \emph{variant generation} (Section \ref{sec:realtimevariant}), with color coding (Section \ref{sec:colors}). This is ideal when up to a few tens of materials are sought and strict control is required over the output (Figures \ref{fig:teaser} (left) and \ref{fig:variant-gen}).
\end{enumerate}
%\vspace{-5mm}
\Revision{The final scene with the newly assigned material models is then to be rendered offline.} In Section \ref{sec:overview} we first present the three learning algorithms, which can be used for the automatic workflow by themselves, while in Section \ref{sec:comb} we show how to combine them to provide an interactive system. %Table \ref{tab:notation} summarizes the notation used throughout the paper.

\section{Learning Algorithms for Material Synthesis}
\label{sec:overview}
 
%\begin{enumerate}
%\item Assign scores, obtain recommendations and visualize them via neural rendering.
%	\begin{enumerate}
%		\item If the recommendations are adequate - accept and proceed,
%		\item If the recommendations point to the right direction, but need refinement - back to step 1., and assign scores to the proposed gallery,
%		\item If the recommendations are not acceptable - back to step 1., and re-assign the scores.
%	\end{enumerate}
%\item Choose workflow.
%	\begin{enumerate}
%		\item Manual workflow - assign the materials to the scene manually and fine-tune a chosen subset of them via \emph{variant generation}. This ideal when up to a few tens of materials are sought and strict control is required over the outputs (Fig. \ref{fig:variant-gen}).
%		\item Automatic workflow - randomly assign the materials to the scene (Figures \ref{fig:teaser} (right) and \ref{fig:wrecking}). Ideal for mass-scale material synthesis when hundreds of material models are sought. 
%	\end{enumerate}
%\end{enumerate}

In this section, we outline the three main pillars of our system: \textbf{Gaussian Process Regression} to perform material learning
and recommendation, a \textbf{Convolutional Neural Network} variant for real-time image prediction, and the \textbf{Gaussian Process
Latent Variable Model} to embed our high-dimensional shader inputs into a 2D latent space to enable the fine-tuning of a select set
of recommended materials (Fig. \ref{fig:architecture}). Table \ref{tab:notation} summarizes the notation used throughout the paper.

\subsection{Material Learning and Recommendation}
\label{sec:gpr}

%\begin{figure*}[ht!]
%\includegraphics[width=1.0\textwidth]{figures/gpr-jsd.png}
%\includegraphics[width=0.35\textwidth]{figures/plots/jsd-time.png}
%\caption{Approximating a simple 1D function with Gaussian Process Regression using 82 na\"ive samples (left). The quality of
the reconstruction is greater than most other regression algorithms (e.g., neural networks), however, as the execution time
scales unfavorably with the number of samples $(O(N^3))$, it should be used in scenarios where samples are scarce. This makes
GPR a suitable choice to learn the material model galleries provided by the user.}
%\label{fig:gpr}
%\end{figure*}

%In this section, we describe an application of Gaussian Process Regression for learning material spaces and show why classical
techniques fail to accurately solve this problem, followed by a way to easily exert control on the output recommendations.
In this section, we propose a combination of Gaussian Process Regression, Automatic Relevance Determination and Resilient
Backpropagation to efficiently perform material learning and recommendation. We also show that these recommendations are easily
controllable.

\paragraph{Material learning.} Gaussian Process Regression (GPR) is a kernel-based Bayesian regression technique that leverages
prior knowledge to perform high-quality regression from a low number of samples\footnote{Throughout this manuscript, we will
use the terms \emph{samples} and \emph{observations} interchangeably.}. 
%As it is endowed with additional information in the form of a covariance kernel, the quality of the reconstruction is greater
than most regression techniques (e.g., neural networks, random forests or SVRs) at low sample counts. 
%However, as several steps contain computing the inverse of the covariance matrix $\covmx$ (Equations (\ref{eq:eq_start})-(\ref{eq:eq_end}))
the execution time scales unfavorably with the number of observations, i.e., $O(n^3)$, it should be used in scenarios where
samples are scarce. This makes GPR a suitable choice to minimize user interaction and learn the material model galleries provided
by the user.
%The goal is to approximate a preference function $\preffungprofparam$ with $\preffunofparam$ from a discrete set of $n$ observations
It can be used to approximate a preference function $\preffunofparam$ from a discrete set of $n$ observations
$ \preffunvec = 
\begin{bmatrix}
    u(\myvec{\param}_1), u(\myvec{\param}_2), \dots, u(\myvec{\param}_n) \\
\end{bmatrix}^T
$,
each of which can be imagined as point samples of a Gaussian where $\myvec{\param}_i \in \mathds{R}^m$ encode the parameters that yield a BSDF model. 
%In our case, $m=19$ \Revision{refers to our parameter space that spans the most commonly used materials, i.e., a combination of diffuse, specular, glossy, transparent and translucent materials. This description is similar to Disney's ``principled shader'' \cite{burley2012physically} to accommodate procedurally textured albedos and displacements (see Section \ref{sec:results} and the supplementary materials for this extended version where $m=38$)}. The goal is to produce $\preffunofparamunknown$ for an arbitrary unknown test input $\paramunknown$ by defining a Gaussian Process over $\preffunvec$.
\Revision{We created a parameter space similar to Disney's ``principled shader'' \cite{burley2012physically} that comes in two versions: the $m=19$ variant spans the most commonly used materials, i.e., a combination of diffuse, specular, glossy, transparent and translucent materials where the extended $m=38$ version additionally supports procedurally textured albedos and displacements (see Section \ref{sec:results} and the supplementary materials).} A Gaussian Process is given by its mean and covariance $\covfun(\param,\param')$ that describes the relation between individual material samples $\param$ and $\param'$. A squared exponential covariance function is given as
\begin{equation}
	\covfun(\param,\param') = \sigma_f^2 \, \exp \Bigg[ -\frac{(\param - \param')^2}{2l^2} \Bigg] + \beta^{-1}\delta_{xx'},
	\label{eq:sqexp}
\end{equation}
with a given $\sigma_f^2$ variance and length scale $l$ where $\beta^{-1}\delta_{xx'}$ is an additional noise term enabled by the Kronecker delta to yield a positive definite covariance matrix. This means that a highly correlated pair $\left\lVert\param - \param' \right\rVert \approx 0$ yields the maximum of the function, i.e., $\covfun(\param,\param') \approx \sigma_f^2 + \beta^{-1}\delta_{xx'}$, leading to a smooth function approximation as $\preffunofparam \approx \preffun(\myvec{\param}')$.
Conversely, if $\left\lVert\param - \param' \right\rVert$ is large, the two observations show negligible correlation, therefore for a rapidly decaying covariance function, $\covfun(\param, \param') \approx 0$.

\begin{figure}[t]
\includegraphics[width=0.475\textwidth]{figures/architecture.png}
\caption{A high-level overview of our pipeline: GPR is used to learn the user-specified material preferences and recommend new materials which are subsequently visualized using our Convolutional Neural Network. Optionally, GPLVM can be used to provide an intuitive 2D space for variant generation.}
\label{fig:architecture}
\end{figure}

The covariance matrix $\covmx$ is given by all possible combinations of the point samples
\begin{equation}
\covmx =
	\begin{bmatrix}
		\covfun(\param_1, \param_1) & \covfun(\param_1, \param_2) & \ldots & \covfun(\param_1, \param_n) \\
		\covfun(\param_2, \param_1) & \covfun(\param_2, \param_2) & \ldots & \covfun(\param_2, \param_n) \\
		\vdots & \vdots & \ddots & \vdots \\
		\covfun(\param_n, \param_1) & \covfun(\param_n, \param_2) & \ldots & \covfun(\param_n, \param_n)
	\end{bmatrix},
\end{equation}
the diagonal of $\covmx$ is therefore always $\sigma_f^2 + \beta^{-1}\delta_{xx'}$. The covariances for the unknown sample $\paramunknown$ are written as 
\begin{equation}
\kstar = 
\begin{bmatrix}
    \covfun(\paramunknown, x_1), \covfun(\paramunknown, x_2), \dots, \covfun(\paramunknown, x_n)
\end{bmatrix}^T,
\end{equation}
(where $x_i \!\in\! \paramunknown$) and $\kstars = \covfun(\paramunknown, \paramunknown)$. We define a zero-mean Gaussian Process over 
$
\begin{bmatrix}
 \preffunvec, \preffun(\paramunknown)
\end{bmatrix}^T
$
with the covariance function $\covfun(\param, \param')$:

\begin{equation}
	\begin{bmatrix}
	 \preffunvec \\ 
	 \preffun(\paramunknown)
	\end{bmatrix}
	\sim \mathcal{N} \Bigg( 0, 
	\begin{bmatrix}
	 \covmx & \kstart \\ \kstar & \kstars
	\end{bmatrix}
	 \Bigg).
\end{equation}

\begin{table}[tb]
\caption{Notation used throughout this paper, in order of occurrence.}
\renewcommand{\arraystretch}{1.075}
\centering
	\begin{tabular}{crr}
		  \toprule
		  Symbol & Description & Type \\
		  \midrule
		  $\myvec{x}$ & BSDF description & Vector \\
		  $\preffungprofparam$ & Preference function (Ground truth) & Scalar \\
		  $\preffunofparam$ & Preference function (GPR prediction) & Scalar \\
		  $n$ & Number of GPR samples & Scalar \\
		  $\paramunknown$ & Unknown BSDF test input & Vector  \\
		  $\preffunvec$ & GPR training set & Matrix \\
		  $m$ & Input BSDF dimensionality & Scalar \\
		  $\covfun(\param,\param')$ & Covariance function & Scalar \\
		  $\sigma_f^2$ & Variance & Scalar \\
		  $l$& Length scale & Scalar \\
		  $\beta^{-1}$& Noise term & Scalar \\
		  $\delta_{xx'}$& Kronecker delta & Scalar \\
		  $\covmx$& GPR covariance matrix & Matrix \\
		  $\boldsymbol{\theta}$& Covariance function parameterization & Vector \\
		  $\bm{\alpha}$& Learning rate & Vector \\
		  $\phi(\myvec{x}^{*})$& CNN image prediction of a BSDF & Matrix \\
		  $\myvec{X}$ & GPLVM training set & Matrix\\
		  $\myvec{L}$ & Low dimensional latent descriptor & Matrix \\
		  $l$ & Latent space dimensionality & Scalar \\
		  $\mathcal{N}( x | \, \mu, \sigma^2 )$ & $\frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(x-\mu)^2}{2\sigma^2})$ & Scalar \\
		  $z$& Number of GPLVM samples & Scalar \\
		  $\myvec{K}^{'}$& GPLVM covariance matrix & Matrix \\
		  $\psi(\myvec{l}^{*})$& Mapping from latent to observed space & Vector  \\
		  $m(\param)$ & $\frac{1}{2}(\preffunofparam + \preffungprofparam)$ & Scalar \\
		  $\tau$ & Recommendation threshold & Scalar \\
		  $r$ & Grid resolution & Scalar  \\
		  $s(\myvec{x}^{*}, \myvec{x}')$ & Similarity (CNN prediction) & Scalar \\
		  $u(\myvec{x}')$& Preference score (GPR prediction) & Scalar \\
		  \bottomrule
	\end{tabular}
\label{tab:notation}
\end{table}

We seek $\preffunofparamunknown$ leaning on the knowledge that the conditional probability $P(\preffunofparamunknown \, | \, \preffunvec )$ follows a Gaussian distribution, therefore the closed-form solution for $\preffunofparamunknown$ and its variance is obtained by
\begin{align*}
	\preffunofparamunknown &= \kstart \covmx^{-1} \preffunvec, \\
	\sigma(\preffunofparamunknown) &= \kstars - \kstar \covmx^{-1} \kstart. \numberthis
	\label{eq:eq_start}
\end{align*}
The quality of the regression depends on the choice and the parameterization of the covariance function. If $\boldsymbol{\theta} = \{\sigma_f^2,l\}$ in (\ref{eq:sqexp}) is chosen poorly, the result will suffer from severe over- or underfitting. To avoid this, a model selection step is performed to maximize the log-likelihood of the observed samples by choosing the appropriate hyperparameters, i.e.,
\begin{equation}
	\log P\big( \preffunvec | \param, \boldsymbol{\theta}\big) = - \frac{1}{2} \preffunvec^T \covmx^{-1} \preffunvec - \frac{1}{2} \log |\covmx| - \frac{n}{2} \log 2\pi.
\end{equation}

The derivatives of the log-likelihood with respect to the hyperparameters and a learning rate $\bm{\alpha}$ are given by
\begin{equation}
\frac{\partial}{\partial \theta_j} \log P\Big( \preffunvec | \param, \boldsymbol{\theta}\Big) = - \frac{1}{2} \tr \Bigg\{ \left( \bm{\alpha} \bm{\alpha}^T \covmx^{-1} \right) \frac{\partial \covmx}{\partial \theta_j} \Bigg\},
	\label{eq:eq_end}
\end{equation}
which can be used with gradient-based optimization techniques. In some cases, the user is looking for a class of material models where, for instance, the surface reflectance properties are of great importance and the choice of albedos is irrelevant (e.g., carpaint material variants). In this case, using one fixed length scale for all features leads to poor predictions. To this end, we have also used Automatic Relevance Determination \cite{mackay1996bayesian,neal2012bayesian} and assigned a $\boldsymbol{\theta}$ to each dimension (with appropriate modifications to (\ref{eq:sqexp})) to open up the possibility of discarding potentially irrelevant features. However, this substantially increases the dimensionality of the optimization problem, to a point where classical methods such as L-BFGS-B \cite{byrd1995limited} and the Scaled Conjugate Gradient method \cite{moller1993scaled} prove to be ineffective: lower learning rates are theoretically able to find the desired minima, but slow down considerably in shallow regions while larger learning rates introduce oscillation near the minima.
To this end, instead of using a fixed step size that is proportional to the local gradient, we adapt the learning rate based on the topology of the error function by using Resilient Backpropagation  \cite{riedmiller1992rprop,blum2013optimization}, a technique originally designed for training neural networks. This significantly reduces the number of required learning steps and has always succeeded finding usable minima in our experiments.

\paragraph{Material recommendation.}
Given enough learning samples, $\preffun$ will resemble the true user preferences, therefore high-scoring material recommendations can be obtained by simply rejection-sampling it against a minimum acceptable score threshold. The rejection rates depend on the properties of $\preffun$: choosing a high threshold will result in high-quality recommendations at the cost of higher rejection rates and decreased variety. Conversely, a larger variety of recommendations can be enforced by lowering the target threshold. Upon encountering a rejection ratio that is unacceptably high, we extend the rejection sampler with a low number of stochastic hillclimbing steps. This offers a considerably improved ratio at the cost of a minimal increase in sample correlation. Generally, we have found setting the recommendation threshold between $40$\% and $70$\% of the maximum possible score to be a good tradeoff between sample quality and variety (Fig. \ref{fig:wrecking}, higher quality (left) -- $70$\%, more variety (right) -- $40$\%).

\begin{figure}[t]
\includegraphics[width=0.475\textwidth]{figures/cnn-arch.png}
\caption{Our decoder network takes the shader description as an input and predicts its appearance. Due to the fact that we have an atypical problem where the input shader dimensionality is orders of magnitude smaller than the output, the input signal is subjected to a series of 1D convolution and upsampling steps.}
\label{fig:cnn_architecture}
\end{figure}

\subsection{Neural Networks and Rendering}
\label{sec:cnn}

%Every image in the recommendation gallery has to be rendered before being presented to the user. This is done via a global illumination, where delaying the artist's work by 40-60 seconds per image for a larger workflow is clearly infeasible. 
After the learning and recommendation step have taken place, a gallery is generated with a prescribed amount of recommendations. The GPR and the recommendation steps take only a few seconds (Table \ref{tab:timings}), however, rendering all 300 recommendations (a typical number throughout our experiments) would take over 4 hours, which is a prohibitively long time for practical use. To cut down the time between the two steps, we introduce a neural network-based solution. Deep neural networks are universal function approximators that have proven to be remarkably useful for classification and regression problems. Typically, the dimensionality of the input is orders of magnitude larger than that of the output (e.g., image classification). Convolutional Neural Networks \cite{lecun1998gradient} (CNNs) excel at solving problems of this form; their advantages include augmenting neurons with a receptive field to take advantage of the locality of information in images and their pooling operations that reduce the number of parameters in each layer.
In this work, we are concerned with an atypical adjoint problem where images are to be predicted pixel by pixel from the shader input, $\phi\!: \mathds{R}^m \!\! \rightarrow \! \mathds{R}^p$, where the input dimensionality $m$ is in the order of tens, which is to be mapped to an output of $p$ dimensions, which represents the largest possible output that fits into the GPU memory along with its decoder network, i.e., in our case, a $410^2$ image with three color channels. We will refer to problems of this kind as \emph{neural rendering}.

\begin{figure*}[ht]
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/best1_captioned.png}
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/worst1_captioned.png}
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/best2.png}
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/worst2.png}
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/best3.png}
\includegraphics[width=0.495\textwidth]{figures/nn-predictions/worst3.png}
\caption{The best (left side) and worst-case (right side) predictions by our neural network on a set of 250 images. Mean PSNR: 37.96dB, minimum: 26.05dB, maximum: 48.70dB.}
\label{fig:nn_predictions}
\end{figure*}

Neural rendering remains an unsolved problem in general. However, as our case is constrained to material modeling, the geometry and lighting setups can be kept constant, therefore it is possible to create a sufficiently deep network and training set to predict images that are nearly indistinguishable from the ones rendered until convergence with full global illumination. Using deconvolutional layers \cite{zeiler2010deconvolutional,noh2015learning} would be a natural choice for $\phi$, however, as few of the most common software packages support it and unwanted artifacts may appear during image generation \cite{odena2016deconvolution}, instead, the input signal is subjected to a series of 1D convolutions, each followed by an upsampling layer to inflate the number of parameters as we advance deeper into the network (Fig. \ref{fig:cnn_architecture}).

\sloppy
This architecture is similar to the decoder part of Convolutional Autoencoders \cite{masci2011stacked} and can be described with the shorthand notation of $4$x$\{$Conv1D($64,3,1$) -- Upsampling($2$)$\}$ -- FC($1000$) -- FC($410^2\!\cdot\!3$), where the parameters of the convolutional layer are \emph{number of filters}, \emph{spatial kernel size} and \emph{strides}, respectively. These layers use exponential linear units \cite{clevert2015fast} with Glorot-initialization \cite{glorot2010understanding} and are trained via the Adam optimizer \cite{kingma2014adam} ($lr\!\!=\!\!10^{-3}, \beta_1\!\!=\!\!0.9, \beta_2\!\!=\!\!0.999, \epsilon\!\!=\!\!10^{-8}, \text{decay}\!\!=\!\!0$). Normally, a network of this size introduces severe overfitting, even in the presence of $\mathcal{L}_1/\mathcal{L}_2$ regularization \cite{nowlan1992simplifying,zou2005regularization} or dropout \cite{srivastava2014dropout}, especially if learning takes place on a given, fixed dataset. However, as we can create our own dataset, i.e., a potentially infinite number of shader-image pairs via rendering, we are able to evade this problem by creating a sufficiently large training set.
%These images were tone mapped using a Gaussian kernel ($\sigma\!=\!1.5$)} 
\Revision{In the interest of efficiency, we have generated $45000$ LDR shader-image pairs with a spatial resolution of $410^2$ and $250$ samples per pixel over 4 weeks on a consumer system with a NVIDIA GeForce GTX TITAN X GPU}, and since no means were required to prevent overfitting, our training process converged to \Revision{$10^{-2}$ (RMSE)}, a negligible but non-zero $\mathcal{L}_2$ training and validation loss (where a zero loss would mean containing all the noise from the training set) in less than 30 hours. \Revision{Our validation set contained $2500$ images and the measured losses correlated well with real-world performance.} This pre-trained network can be reused as long as the shader description remains unchanged and the inference of a new image typically takes $3$ to $4$ milliseconds. A further advantage of this architecture is that similarly to Denoising Autoencoders \cite{vincent2010stacked}, it also performs \emph{denoising} on the training samples.

A key observation is that the more layers the neural network contains, the more high-frequency details it will be able to represent (a similar effect has been observed in Fig 4., Saito et al. \shortcite{saito2016photorealistic}). This means that our proposed neural network contains enough layers to capture the important features to maximize visual quality, but not enough to learn the high-frequency noise contained in the dataset. 
\setlength{\columnsep}{6.5pt}
\begin{wrapfigure}{hr}{0.2\textwidth}
  \centering
  %\hspace{-2pt}
  \vspace{-12pt} % With some placements, this is useful
  %\vspace{+10pt} % 
  \includegraphics[width=0.2\textwidth]{figures/noise_filter.png}
  %\vspace{-20pt} % With some placements, this is useful
  %\vspace{-23.75pt} % 
  \vspace{-26.5pt} % 
%  \vspace{-28.5pt} %
%  \caption{Caption}
\end{wrapfigure}
%This is a remarkable scenario where the network's predictions are closer to the ground truth images than the training samples used to train it. 
Normally, this denoising process requires the presence of auxiliary feature buffers and longer computation times \cite{sen2012filtering} or other deep learning approaches using more complex architectures \cite{bako2017kernel}. In our case, this step does not require any additional complexity and is inherent to the structure of our network. Beyond producing less noisy images in real time, this is also a significant advantage in easing the computational load of creating new datasets as the training images do not have to be rendered until convergence. We took advantage of this by using only $250$ samples per pixel for each training image, which took six times less than the standard $1500$ samples that would be required for perfectly converged training samples, cutting down the total time to produce the training set from $24$ to $4$ weeks.
%\begin{figure}[ht!]
%\includegraphics[width=0.475\textwidth]{figures/noise_filter.png}
%\caption{Caption.}
%\label{fig:nn-noise}
%\end{figure}

\subsection{Latent Space Variant Generation}
\label{sec:gplvm}
After being presented with a gallery of material models, the user may find that some recommendations are close to their preference, but some require fine-tuning to fit their artistic vision for a scene. Adjusting the materials by hand requires domain expertise and a non-trivial amount of trial and error. In our solution, we seek a dimensionality reduction technique that maps the shader inputs into a 2D latent space where similar materials can be intuitively explored. Non-linear manifold learning techniques come as a natural choice for this kind of problem, however, the most well-known methods \cite{belkin2003laplacian,maaten2008visualizing,tenenbaum2000global} could not find coherent structures in 2D. Beyond that, these techniques are unable to interpolate between the embedded samples, therefore they would only be useful for visualization in a latent space, but not for exploration. Dimensionality reduction with autoencoders \cite{hinton2006reducing} would require orders of magnitude more samples to work properly (we have at most a few tens at our disposal), which is infeasible as it would be too laborious for the user to provide that many scores through the presented gallery.

The Gaussian Process Latent Variable Model (GPLVM) \cite{lawrence2004gaussian} is a non-linear dimensionality reduction technique which is able to embed a few tens of high-scoring materials from the gallery
$ \myvec{X} = 
\begin{bmatrix}
 \cdots \myvec{x}_i \cdots \\
\end{bmatrix}^T
$
with $\myvec{x}_i \in \mathds{R}^m$ into a set of low dimensional latent descriptors
$ \myvec{L} = 
\begin{bmatrix}
 \cdots \myvec{l}_i \cdots \\
\end{bmatrix}^T
$ with $\myvec{l}_i \in \mathds{R}^l$.
Typically, $m \! \gg \! l$, in our case, $m=19$ and $l=2$ to make sure that variant generation can take place conveniently on a 2D plane. The likelihood of the high-dimensional data using $z$ high-scoring training samples is given as
\begin{equation}
%P \Big( \myvec{X} | \myvec{L}, \boldsymbol{\theta} \Big) = \prod_{z=1}^{Z} \mathcal{N} \Big( \myvec{X}_z | \, 0, \; \myvec{K}(\myvec{L},\myvec{L}|\boldsymbol{\theta}) + \sigma^2 \myvec{I} \Big).
P \Big( \myvec{X} | \myvec{L}, \boldsymbol{\theta} \Big) = \prod_{i=1}^{z} \mathcal{N} \Big( \myvec{X}_i | \, \myvec{0}, \; \myvec{K}^{'} \!\! + \beta^{-1} \myvec{I} \Big),
\end{equation}
where $\myvec{K}^{'}$ is a covariance matrix similar to $\covmx$ containing elements akin to (\ref{eq:sqexp}) with a substitution of $\covfun(\param,\param') \rightarrow k(\myvec{l},\myvec{l}^{'})$. In this case, the optimization takes place jointly over the latent values in $\myvec{L}$ and $\boldsymbol{\theta}$, i.e.,
\begin{equation}
	\myvec{L}^*, \boldsymbol{\theta}^* = \argmax_{\myvec{L},\boldsymbol{\theta}} \, \log \Big[ P \Big( \myvec{X} | \myvec{L}, \boldsymbol{\theta} \Big) \Big].
\end{equation}
Even though PCA disregards the non-linear behavior of BSDF parameters \cite{lafortune1997non}, it serves as a formidable initial guess for $\myvec{L}$ \cite{lawrence2004gaussian} and $\boldsymbol{\theta}$ is initialized with a wide prior. Beyond the ability to learn efficiently from a handful a samples, a further advantage of this method is that a new mapping can be made between the latent and observed space $\myvec{x}^{*}=\psi(\myvec{l}^{*})$, i.e.,
\begin{equation}
	\begin{bmatrix}
	 \myvec{X} \\ 
	 \myvec{x}^{*}
	\end{bmatrix}
	\sim \mathcal{N} \Bigg( 0, 
	\begin{bmatrix}
	 \myvec{K}^{'} & \myvec{k}_{*}^{'T} \\ \myvec{k}_{*}^{'} & \myvec{k}_{**}^{'}
	\end{bmatrix}
	 \Bigg),
\end{equation}
yielding the final closed-form solution
\begin{align*}
	\psi(\myvec{l}^{*})  &= \myvec{k}_{*}^{'T} \myvec{K}^{'-1} \myvec{X}, \\
	\sigma(\psi(\myvec{l}^{*})) &= \myvec{k}_{**}^{'} - \myvec{k}_{*}^{'} \myvec{K}^{'-1} \myvec{k}_{*}^{'T}. \numberthis
	\label{eq:eq_gplvm}
\end{align*}
By storing and reusing $\myvec{K}^{'-1}$, this mapping can be done in negligible time, which allows the user to rapidly generate new material variants in this 2D latent space.
