# Deep Learning Partial Least Squares

Nicholas Polson111Email<EMAIL_ADDRESS>
Booth School of Business
University of Chicago
Vadim Sokolov
Department of Systems Engineering
and Operations Research
George Mason University
Jianeng Xu
Booth School of Business
University of Chicago

###### Abstract

High dimensional data reduction techniques are provided by using partial least
squares within deep learning. Our framework provides a nonlinear extension of
PLS together with a disciplined approach to feature selection and architecture
design in deep learning. This leads to a statistical interpretation of deep
learning that is tailor made for predictive problems. We can use the tools of
PLS, such as scree-plot, bi-plot to provide model diagnostics. Posterior
predictive uncertainty is available using MCMC methods at the last layer. Thus
we achieve the best of both worlds: scalability and fast predictive rule
construction together with uncertainty quantification. Our key construct is to
employ deep learning within PLS by predicting the output scores as a deep
learner of the input scores. As with PLS our X-scores are constructed using
SVD and applied to both regression and classification problems and are fast
and scalable. Following Frank and Friedman (1993), we provide a Bayesian
shrinkage interpretation of our nonlinear predictor. We introduce a variety of
new partial least squares models: PLS-ReLU, PLS-Autoencoder, PLS-Trees and
PLS-GP. To illustrate our methodology, we use simulated examples and the
analysis of preferences of orange juice and predicting wine quality as a
function of input characteristics. We also illustrate Brillinger’s estimation
procedure to provide the feature selection and data dimension reduction.
Finally, we conclude with directions for future research.

Key Words: Deep Learning, Partial Least Squares, PLS-ReLU, PLS-Autoencoder,
PLS-Trees, PLS-GP, Principal Component Analysis, Nonlinear Classification,
Trees, Gaussian Processes, Bayesian learning, Shrinkage, Dimension Reduction

##  1 Introduction

Many high dimensional regression problems involve identifying a functional
that maps a high-dimensional input variable $x\in R^{p}$ to an output variable
$y\in R^{q}$,

$y=\phi(x)+\epsilon.$

The function $\phi$ is estimated using the observed matrices of inputs $X\in
R^{n\times p}$ and corresponding outputs $Y\in R^{n\times q}$. A high-
dimensional input $x$ naturally arise in many applications, such as marketing,
genetics and robotics. For example, in business applications, it is typical to
increase the dimensionality of the input $x$ by introducing new hand-coded
predictors (Hoadley, 2001). This expanded input can include terms such as
interactions, dummy variables or nonlinear functionals of the original
predictors.

Many high dimensional statistical models require data dimension reduction
methods. Following Breiman (2001) we represent data as being generated by a
black box with an input vector $x$. The black box generates output $y$, or a
predictive distribution over the output $p(y\mid x)$ that describes
uncertainty in predicting $y$ from $x$. Fisher (1922) and Cook (2007) describe
clearly the issue of dimension reduction. Although it was typical to find
predictors via screening and plotting them against the output variable Li and
Nachtsheim (2007), Fisher argued with this approach and said in a classic
empirical analysis, “The meteorological variables to be employed must be
chosen without reference to the actual crop record”. Many researchers (Cox,
1968) since then showed counterexamples to Fisher’s approach, when
dimensionality reduction techniques that do not use the output variable lead
to poor results Fearn (1983). In particular, the output variables should be
used for dimension reduction when the input variables are controlled by the
experimenter (Cox and Reid, 2000).

Partial least squares (Wold (1975) and Wold et al. (1984)) has been applied in
many fields, in particular, chemo-kinetics, but has received little attention
in machine learning. On the other hand, nonlinear dimensionality reduction
techniques for feature selection such as deep learning are prevalent. Our goal
here is to incorporate deep learning into partial least squares to provide a
framework for statistical interpretation and predictive uncertainty in high
dimensional problems. Thus we provide a very flexible nonlinear extension of
PLS that merges the two cultures of black box and statistical modeling as
proposed by Breiman (2001). In doing so, we introduce a variety of new partial
least squares models: PLS-ReLU, PLS-Autoencoder, PLS-Trees and PLS-GP.

Our work builds on the literature of nonlinear PLS (Wold et al., 1984; Frank,
1994; Malthouse et al., 1997; Eriksson et al., 2009; Liu et al., 2019). Deep
Learners model $\phi$ as a superposition of univariate affine functions
(Polson and Sokolov, 2017) and are universal approximators. Whilst Gaussian
Process (Gramacy and Lee, 2008; Higdon et al., 2008) are also universal
approximators and can capture relations of high complexity, they typically
fail to work in high dimensional settings. Tree methods can be very effective
in high dimensional problems. Hierarchical models are flexible stochastic
models but require high-dimensional integration and MCMC simulation. Deep
learning, on the other hand, is based on scalable fast gradient learning
algorithms such as Stochastic Gradient Descent (SGD) and its variants. Modern
computational techniques such as automated differentiation (AD) and
accelerated linear algebra (XLA) are available to perform stochastic gradient
descent (SGD) at scale (Polson and Sokolov, 2020), thus avoids the curse of
dimensionality by simply pattern matching and using interpolation to predict
in other regions of the input space. User friendly libraries such as
TensorFlow or PyTorch provide high level interfaces and make deep learning
modeling accessible without requiring to implement any of the low level linear
algebra routines. The algorithmic culture has achieved much success in high
dimensional problems.

Rather than using shallow additive architectures as in many statistical
models, deep learning uses layers of semi affine input transformations to
provide a predictive rule. Applying these layers of transformations leads to a
set of attributes (a.k.a features) to which predictive statistical methods can
be applied. Thus we achieve the best of both worlds: scalability and fast
predictive rule construction together with uncertainty quantification. Sparse
regularisation with un-supervised or supervised learning finds the features.
We clarify the duality between shallow and wide models such as PCA, PPR and
deep but skinny architectures used in deep learning such as auto-encoders,
MLPs, CNN, and LSTM. The connection with data transformations is of practical
importance for finding good network architectures.

There are a number of advantages of our methodology over traditional deep
learning. First, PLS uses linear methods to extract features and weight
matrices. Due to a fundamental estimation result of Brillinger (2012), we show
that these estimates can provide a stacking block for our DL-PLS model.
Moreover, traditional diagnostics (bi-plots, scree-plots) to diagnose the
depth and nature of the activation functions allow us to find a parsimonious
model architecture. On the theoretical side, we follow Frank and Friedman
(1993) and provide a Bayesian shrinkage interpretation of our model. This
allows us to discriminate between unsupervised and supervised learning method.

On the empirical side, we illustrate our methodology using simulated examples
and analysis of preferences of orange juice data and predicting wine quality
as a function of input variables. We also illustrate the novel estimation
procedure suggested by Brillinger to provide the feature selection and data
dimension reduction. Lindley (1968) provides an alternative fully Bayesian
decision theoretic approach to select the number of components in the model.
Hahn et al. (2013) provide a full Bayes analysis of stochastic factor models.

The rest of the paper is outlined as follows. Section 1.1 provides
definitions, notation and connections with previous work. Section 2 describes
partial least squares from the perspective of Bayesian shrinkage. Section 3
provides our general modeling strategy of deep learning within partial least
squares. We provide a nonlinear extension of partial least squares by using
deep layer to model the relationship between the scores generated in PLS.
Section 4 provides applications, We illustrate a key theoretical result of
Brillinger on how to estimate deep layered models and show how to recover ReLU
and tanh activation functions. The orange juice dataset is re-analysed and we
find evidence for nonlinearity in the scores architecture. Section 5 concludes
with directions for future research.

###  1.1 Connection with Previous Work

From a probabilistic view point, it is natural to view input-output pairs as
being generated from some joint distribution

$(y_{i},x_{i})\sim p(y,x),\;i=1,\dots,N.$

The joint distribution $p(y,x)$ is characterised by its two conditionals
distributions

  * •

$p(y\mid x)$, which can be used for prediction $\hat{y}=E(y\mid x)$ together
with full predictive uncertainty quantification using the probabilistic
structure.

  * •

$p(x\mid y)$. This conditional is high-dimensional and we need to perform
dimension reduction and find an efficient set of nonlinear features to be used
as predictors in step 1. Selection using sparsity and deep learners are
typical methods used.

In our PLS-DL model we use a hierarchical linear regression model as the last
layer. This allows us to perform a full Bayes inference at the level, namely
calculating the predictive $p(y\mid x)$ helps with assessing prediction
uncertainty. Using stochastic regressors allows us to use inference methods
proposed by Brillinger (2012).

There are several types of projection-based dimensionality deduction
techniques in regression. All of those rely on finding orthogonal projections
of columns of $X$ and (possibly) $Y$ into a lower dimensional liner space.
Essentially, the problem us to find a new basis and to find the projections
into this basis. The basis vectors are called loadings and the projection is
called the score.

Principal Components is arguable the most widely used dimensionality reduction
technique in high-dimensional settings. The dimensionality reduction is
achieved by changing a basis and using singular vectors of the input matrix
$X$ as new basis vectors. The singular vector decomposition $X=USV^{T}$, where
$V$ is $p\times p$ orthogonal matrix, $U$ is $n\times n$ orthogonal and $S$ is
$n\times p$ diagonal matrix, provides the new set of basis vectors
(directions) $V$ and the projections into this new basis $XV=US$ are called
the principal components (scores). Then dimensionality reduction is achieved
by projecting inputs into space spanned by the first $L<p$ right singular
vectors $V_{L}=(v_{1},\ldots,v_{L})$ which correspond to the $L$ largest
singular values $s_{1}\leq s_{2},\ldots,\leq s_{L}$. The resulting projection
leads to the best $L$-rank approximation of the input matrix $X$, $X\approx
U_{L}S_{L}V_{L}^{T}$ (Golub and Van Loan, 2013).

Principal component dimension reduction is independent of $y$ and can easily
discard information that is valuable for predicting the desired output. Sliced
inverse regression (SIR) (Li, 1991; Cook, 2009; Jiang and Liu, 2014) overcomes
this drawback and assumes that conditional distribution can be indexed by $L$
linear combinations

$p(y\mid x)=p(y\mid Wx),$

and by computing the set of $L$ basis vectors (rows of $W$) using both $x$ and
$y$.

Partial Least Squares (PLS) also uses both $X$ and $Y$ to calculate the
projection, further it simultaniously finds projections for both input $x$ and
output $y$, it make it applicable to the problems with high-dimensional output
vector as well as input vector. Let $Y$ be $n\times q$ matrix of observed
outputs and $X$ be an $n\times p$ input matrix. PLS finds a projection
directions that maximize covariance between $X$ and $Y$, the resulting
projections $U$ and $T$ for $Y$ and $X$, respectively are called score
matrices, and the projection matrices $P$ and $Q$ are called loadings. The
$X$-score matrix $T$, has $L$ columns, one for each “feature” and $L$ is
chosen via cross-validation. The key principle is that $T$ is a good predictor
of $U$, the $Y$-scores. This relations are summarized by the equations below

$\displaystyle Y$ $\displaystyle=UQ+E$ $\displaystyle X$ $\displaystyle=TP+F$

Here $Q$ and $P$ are orthogonal projeciton (loading) matrices and $T$ and $U$
are $n\times L$ are projections of $X$ and $Y$ respectively.

Originally PLS was developed to deal with the problem of collinearity in
observed inputs. Although, principal component regression also addresses the
problem of collinearity, it is often not clear which components to choose. The
components that correspond to the larges singular values (explain the most
variance in $X$) are not necessarily the best ones in the predictive settings.
Also ridge regression addresses this problem and was criticized by (Fearn,
1983). Further ridge regression does not naturally provide projected
representations of inputs and outputs that would make it possible to combine
it with other models as we propose in this paper. Thus, PLS seem to be the
right method for high-dimensional problems when the goal is to model non-
linear relations using another model.

In the literature, there are two types of algorithms for finding the
projections (Manne, 1987). The original one proposed by Wold et al. (1984)
which uses conjugate-gradient method (Golub and Van Loan, 2013) to invert
matrices. The first PLS projection $p$ and $q$ is found by maximizing the
covariance between the $X$ and $Y$ scores

$\displaystyle\operatorname*{maximize\quad}_{p,q}$
$\displaystyle\left(Xp\right)^{T}\left(Yq\right)$ subject to
$\displaystyle||p||=||q||=1.$

Then the corresponding scores are

$t=Xp,\text{ and }u=Yq$

We can see from the definition that the directions (loadings) for $Y$ are the
right singular vectors of $X^{T}Y$ and loadings for $X$ are the left singular
vectors. The next step is to perform regression of $T$ on $U$, namely
$U=T\beta$. The next column of the projection matrix $P$ is found by
calculating the singular vectors of the residual matrices
$(X-tp^{T})^{T}(Y-T\beta q^{T})$.

The final regression problem is solved $Y=UQ=T\beta Q=XP^{T}\beta Q$. Thus the
PLS estimate is

$\beta_{\mathrm{PLS}}=P^{T}\beta Q.$

Helland (1990) showed that PLS estimator can be calculated as

$\beta_{\mathrm{PLS}}=R(R^{T}S_{xx}R)^{-1}R^{T}S_{xy}$

where $R=(S_{xy},S_{xx}S_{xy},\ldots,S_{xx}^{q-1}S_{xy})$,

$S_{xx}=\dfrac{X^{T}(I-{\bf 11^{T}}/n)X}{n-1},$ $S_{xy}=ave(y\bm{x}).$

Further, Helland (1988) proposed an alternative algorithm to calculate the
parameters
For $K=1,\ldots,p$, set $y_{0}=y,\bm{x}_{0}=\bm{x}$. Let $(y,\bm{x})$ be
centered and standardised. Given $V=ave(\bm{x}\bm{x}^{T})$ and
$\bm{s}=ave(y\bm{x})$, for $K=1$ to $p$ do

  1. 1.

$\bm{s}_{k}=V^{K-1}\bm{s}$

  2. 2.

$\hat{y}_{K}=OLS(y\;{\rm on}\;(\bm{s}_{k}^{T}\bm{x})_{k=1}^{K})$

Stone and Brooks (1990) discuss PLS and cross-validation.

#### Non-Linear Models and Ridge Functions

The main building component of a deep learning model is a ridge function,
which finds one dimensional-structure in the data (Pinkus, 2015). Deep
learning is then a composition of several ridge functions to identify non-
linear low-dimensional structure present in high-dimensional data.

A ridge function is a multivariate function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ of the form $f(x)=g(w^{T}x)$ where $g$
is univariate function and $x,w\in\mathbb{R}^{n}$. The non-zero vector $w$ is
called the direction. The name ridge comes from the fact that the function is
constant along directions that are orthogonal to $w$, for a direction $\eta$
such that $w^{T}\eta=0$, we have

$f(x+u)=g(w^{T}(x+u))=g(w^{T}x)=f(x)$

Arguably, this is the simplest multivariate function. Klusowski and Barron
(2016) provided bounds for the bias when noisy data is modeled by smooth ridge
functions.

A generalization of a ridge function is a function of the form

$f(x)=g(Wx),$

where $W$ is a non-zero $d\times n$ matrix of weights and univariate function
$g$ is applied element-wise. Deep learning is nothing but a composition of
such generalized ridge functions.

Ridge functions are used in high-dimensional statistical analysis. For
example, projection pursuit algorithm approximates the input-output relations
using a linear combination of ridge functions (Friedman and Stuetzle, 1981;
Huber, 1985; Jones, 1992)

$\phi(x)=\sum_{i=1}^{r}g_{i}(w_{i}^{T}x),$

where both the directions $w_{i}$ and functions $g_{i}$ are variables and
$w_{i}^{T}x$ are one-dimensional projections of the input vector. Projection
pursuit regression is a dimensionality reduction technique. The vector
$w_{i}^{T}x$ is a projection of the input vector $x$ onto a one-dimensional
space and $g_{i}(w_{i}^{T}x)$ can be though as a feature calculated from data.
Projection pursuit algorithm is an iterative procedure that finds a predictive
function iteratively by starting with $\phi_{0}=0$ and proceeds

$w_{i+1}\in\arg\min_{w}\lVert\phi_{i+1}-y\rVert_{p};\leavevmode\nobreak\
\leavevmode\nobreak\ \phi_{i+1}=\phi_{i}+g_{i+1}(w^{T}x).$

Nonlinear PLS was first considered by Kramer (1991) and Malthouse et al.
(1997) who proposed projecting (non-orthogonally) the input variables into an
$L$-dimensional surface which is parametrized by a neural network. Those
projections, which are equivalent of scores vectors in PLS are then used for
the prediction of the output. Rosipal and Trejo (2001) develops a kernel PLS
algorithm which is a non-linear regression in high dimensions. First, they
find a transformation of the input vector $x\rightarrow\Phi(X)$ and then use
PLS to map $\Phi(X)$ to $y$. The features $\Phi(x)$ are calculated as using a
kernel function that depends on both $x$ and $y$. It can be shown that kernel
regression is similar to a radial basis function network Haykin (1998).

#### Deep learning

DL overcomes many classic drawbacks by jointly estimating non-linear function
that does dimensionality reduction and prediction using information on $Y$ and
$X$ as well as their relationships, and by using $L>2$ layers.

A deep learning is simply a composition of ridge functions. The prediction
rule is embedded into a parameterised deep learner, a composite of univariate
semi-affine functions, denoted by $G_{W}$ where $W=[w^{(1)},\ldots,w^{(L)}]$
represents the weights of each layer of the network. A deep learner takes the
form of a composition of link functions

$G_{W}=G_{1}\circ\dots\circ G_{L}\;{\rm where}\;G_{L}=g_{L}(w_{L}x+b_{L})$

where $g_{L}$ is a univariate link function.

If we choose to use nonlinear layers, we can view a deep learning routine as a
hierarchical nonlinear factor model or, more specifically, as a generalised
linear model (GLM) with recursively defined nonlinear link functions, see
Polson and Sokolov (2017) and Tran et al. (2020) for further discussion.

Diaconis and Shahshahani (1984) use nonlinear functions of linear
combinations. The hidden factors $Z_{i}=W_{i1}X_{1}+\ldots+W_{ip}X_{p}$
represents a data reduction of the output matrix $X$. The model selection
problem is to choose how many hidden units. We will show that PLS scree-plot
can be used as a guide to find $L$. Predictive cross-validation can also be
sued.

We now turn to estimation of single index models with non-Gaussian regressors.

###  1.2 Projection to Single Index Model

Brillinger (2012) considers the single-index model with non-Gaussian
regressors where $(Y,X)$ are stochastic with conditional distribution

$Y\mid X\sim N(g(\alpha+\beta X),\leavevmode\nobreak\ \sigma^{2}).$

Here $\beta X$ is the single features found by data reduction from high
dimensional $X$. Let $\hat{\beta}_{OLS}$ denote the least squares estimator
which solved $X^{T}Y=X^{T}X\beta$. By Stein’s lemma,

$\mathrm{cov}(Y,X)=\beta\mathrm{cov}(g(\alpha+\beta X),\alpha+\beta
X)\mathrm{var}(X)/\mathrm{var}(\alpha+\beta X)$

Then $\hat{\beta}$ is consistent as

$\hat{\beta}_{OLS}=\mathrm{cov}(Y,X)/\mathrm{var}(X)\rightarrow k\beta\;\;{\rm
where}\;\;k=\mathrm{cov}(g(\alpha+\beta X),\alpha+\beta
X)\mathrm{var}(X)/\mathrm{var}(\alpha+\beta X)$

Hence, $\hat{\beta}_{OLS}$ estimator is proportional to $\beta$. We can also
non-parametrically estimate $g(u)$ by plotting
$(\hat{\beta}x_{j},y_{j}),\leavevmode\nobreak\ j=1,\ldots n$ and smoothing
$y_{j}$ values with $\hat{\beta}x_{j}$ near $u$.

Hence, when the $X$s are Gaussians and independent of the error, we have the
relationship $\mathrm{cov}(Y,X)=k\beta\mathrm{var}(X)$. This relationship
follows from the weaker assumption

$E(Y|X)=g(\alpha+\beta X)$

Hence, this approach can be applied to binary classification with
$Prob(Y=1|X)=g(\alpha+\beta X)$ and other models such as survival models.

Naik and Tsai (2000) apply this to partial least squares. Hall and Li (1993)
show that when dimension of $X$ is high then most nonlinear regressions are
still nearly linear.

The following theorem extends Brillinger’s result to deep learning with ReLU
activation.

######  Theorem 1.

Given a deep learning model of the form

$E(Y|X)=G_{1}\circ\dots\circ G_{L}(X)$

where $G_{1}(Z)=B_{1}Z$ and $G_{l}(Z)=\max(B_{l}Z,0)$ for $l=2,3,...,L$.
$B_{L}\in\mathbb{R}^{D\times p}$ and $B_{l}\in\mathbb{R}^{D\times D}$ for
$l=1,2,...,L-1$. $X\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}^{D}$.

Then we can consistently estimate $(B_{1},\ldots,B_{L})$ recursively using
OLS.

###### Proof.

Define $Z_{L}:=X,\;Z_{0}:=Y$ and $Z_{l-1}:=G_{l}(Z_{l})$ for $l=1,2,3,...,L$.
First, we write the deep learning model as

$E(Y|X)=G^{(L)}(B_{L}Z_{L})\;\;{\rm where}\;\;G^{(L)}=G_{1}\circ\dots\circ
G_{L-1}\circ\max(\cdot,0).$

Running OLS of $Y$ on $X$, allows using Brillinger’s result to consistently
estimate $B_{L}$. That is, we get $\hat{B}_{L}\approx k_{L}B_{L}$ where
$k_{L}$ is a diagonal matrix. Then construct

$\hat{Z}_{L-1}=\max(\hat{B}_{L}X,0)\approx\max(k_{L}B_{L}X,0)=k_{L}Z_{L-1}$

and run OLS of $Y$ on $\hat{Z}_{L-1}$. Note that now the corresponding model
becomes

$E(Y|X)=G^{(L-1)}(B_{L-1}Z_{L-1})=G^{(L-1)}\left(k_{L}^{-1}B_{L-1}\cdot
k_{L}Z_{L-1}\right)\;\;{\rm where}\;\;G^{(L-1)}=G_{1}\circ\dots\circ
G_{L-2}\circ\max(\cdot,0).$

According to Brillinger, the resulting OLS estimate $\hat{B}_{L-1}$ satisfies
$\hat{B}_{L-1}\approx k_{L-1}k_{L}^{-1}B_{L-1}$. And correspondingly,

$\hat{Z}_{L-2}=\max(\hat{B}_{L-1}\hat{Z}_{L-1},0)\approx\max(k_{L-1}B_{L-1}Z_{L-1},0)=k_{L-1}Z_{L-2}.$

Therefore, by induction, we can estimate $(B_{1},\ldots,B_{L})$ recursively
solely using OLS. ∎

Erdogdu et al. (2016) proposes an efficient algorithm for calculating the
proportionality constant $k$. The proportionality was first observed by Fisher
for logistic regression with Gaussian design. Antoniadis et al. (2004) provide
techniques for Bayesian estimation of the single-index models using B-spline
approximate $g$ functions.

Xia et al. (2009) consider the generalized single index model

$Y=g(BX)+\epsilon$

where $B\in\mathbb{R}^{D\times p}$ is orthogonal matrix. $BX$ is known as the
effective dimension reduction (EDR)

Stein’s lemma can also be extended to scale mixtures of Gaussians (Gron et al.
(2012)) and a similar analysis holds.

##  2 Bayesian Shrinkage via Partial Least Squares

###  2.1 Ridge, PCR and PLS

Linear prediction rules without shrinkage, $\hat{Y}=X\beta$ where
$\beta=(X^{T}X)^{-1}X^{T}Y$, have the caveat of high variance due to the
typical ill-conditioned nature of the $(X^{T}X)$-design matrix. Thus shrinkage
methods are been developed.

RR, PCR and PLS are all operationally the same. They bias coefficient vector
away from directions in which the predictors have low sampling variability— or
equivalently, sway from the ”least important” principal components of $X$

  1. 1.

With high collinearity in which variance dominates the bias, the performance
of Ridge, PCR and PLS tend to be quite similar. The three methods shrink
heavily along those directions of small predictor variability. They are also
invariant under rotation, but not scaling.

  2. 2.

Unlike Ridge and PCR, PLS biases the coefficient vector towards directions in
the predictor space that preferentially predict the high spread directions in
the response space. It not only shrinks the OLS in some eigen-directions, but
expands it in others (related to $K$-th eigenvalue).

  3. 3.

For PLS, the scale factors along each of the eigen-directions are not linear
in the response values (smooth but not monotonic). They also depend on the OLS
solution (not the length, but just the relative values).

  4. 4.

PLS typically uses fewer components to achieve the same overall shrinkage as
PCR, generally about half as many components. If the number of component being
equal, PLS has less bias and greater variance, compared with PCR.

  5. 5.

For multi-response problems, two-block PLS or the multi-response ridge doesn’t
dramatically do better than the corresponding uni-response procedures applied
separately to the individual responses.

As with any automatic procedure, there are caveats where certain empirical
datasets can load on the eigen-rvectors with the smallest eigen-values, see,
for example, Fearn (1983), Polson and Scott (2012). PLS can be viewed as an
data pre-processing step before implementing deep learning.

###  2.2 Bayesian Shrinkage Interpretation

We will rely heavily on the Bayesian shrinkage interpretation of PCR and PLS
due to Frank and Friedman (1993). Polson and Scott (2010, 2012) provide a
general theory of global-local shrinkage and, in particular, analyze $g$-prior
and horseshoe shrinkage. PLS behaves differently from standard shrinkage rules
as it can shrink away from the origin for certain eigen-directions.

First, start with the SVD decomposition of the input matrix: $X=UDW^{T}$. If
$n>p$ then of full rank. Here $D=diag(d_{1},\ldots,d_{p})$ nonzero ordered
$e_{1}>\ldots>e_{p}$ singular values. Then $W=(w_{1},\ldots,w_{p})$ is the
matrix of eigenvectors for $S=X^{T}X$. We can then transform the original
first layer to an orthogonal regression, namely $y=(UD)W^{T}\beta$ with
corresponding OLS estimator $\hat{\alpha}=(Z^{T}Z)^{-1}Z^{T}y=D^{-1}U^{T}y$.

#### Ridge, PCR and PLS Algorithm

PCR has a long history in statistics. This is an unsupervised approach to
dimension reduction (no $y$’s). Specifically, we first center and standardise
$(y,\bm{x})$.

Then, we provide an SVD decomposition of

$V=ave(\bm{x}\bm{x}^{T})$

This find the eigen-values $e_{j}^{2}$ and eigenvectors arranged in descending
order, so we can write

$V=\sum_{j=1}^{p}e_{j}^{2}{\bm{v}}_{j}{\bm{v}}_{k}^{T}.$

This leads to a sequence of regression models $(\hat{Y}_{0},...,\hat{Y}_{L})$
with $\hat{Y}_{0}$ being the overall mean and

$\hat{Y}_{L}=\sum_{l=0}^{L}(ave(w_{l}^{T}\bm{x})/e_{l}^{2})\bm{v}_{l}^{T}\bm{x}$

Therefore, PLS finds ”features”
$Z_{K}=\\{\bm{v}_{k}^{T}\bm{x}\\}_{k=0}^{K}=\\{\bm{f}_{k}\\}_{k=0}^{K}$.

The key insight is that all of the estimators are of the from

$\hat{Y}^{M}=\sum_{j=1}^{L}f_{j}^{M}\hat{\alpha}_{j}\bm{v}_{j}^{T}\bm{x}$

where $f_{j}$ are scale factors. $M$ denotes method (e.g. RR, PCR, PLS). $L$
is the rank of $\bm{V}$ (number of nonzero $e_{k}^{2}$). For PCR, the scale
factors are $f_{j}=1$ for top $L$ eigenvectors Frank and Friedman (1993)

$\displaystyle f_{j}^{RR}$
$\displaystyle=e_{j}^{2}/(e_{j}^{2}+\lambda),\text{where $\lambda$ is a fixed
regularization parameter}$ $\displaystyle f_{j}^{PCR}$
$\displaystyle=\begin{cases}1,&e_{j}^{2}\geq e_{L}^{2}\\\
0,&\text{otherwise}\end{cases}$ $\displaystyle f_{j}^{PLS}$
$\displaystyle=\sum_{k=1}^{K}\theta_{k}e_{j}^{2k}\;,\;{\rm
where}\;\theta=w^{-1}\eta,\;\eta_{k}=\sum_{j=1}^{p}\hat{\alpha}_{j}^{2}e_{j}^{2(k+1)}.$

The corresponding shrinkage factors for RR and PCR are typically normalized so
that they give the same overall shrinkage so that the length of solution
vector are the same ($|\hat{\beta}_{RR}|=|\hat{\beta}_{PLS}|$).

This scale factors provide a diagnostic plot: $f_{j}$. If any $f_{j}>1$ then
one can expect supervised learning (a.k.a. PLS with $Y$’s influence the
scaling factors) will lead to different predictions than unsupervised learning
(a.k.a. PCR with solely dependent on $X$). In this sense, PLS is an optimistic
procedure in that the goal is to maximise the explained variability of the
output in sample with the hope of generalising well out-of-sample.

#### Model selection (a.k.a. dimension reduction)

The goal of PCR is to minimize predictive MSE

$\hat{L}=\arg\min_{L}ave(y-\hat{y}_{L})^{2}$

The choice of $L$ is determined via predictive cross-validation. The $l$th
model is simple regression of $y$ on $f_{L}=\bm{v}_{l}\bm{x},l=1,\ldots L$.
Mallows (1973) $C_{p}$ and $C_{L}$ provide the relationship between shrinkage
and model selection.

#### Dropout

This is a model selection technique designed to avoid over-fitting in deep
learning. This is done by removing input dimensions in $X$ randomly with a
given probability $p$. For example, suppose that we wish to minimise MSE,
$\|Y-\hat{Y}\|^{2}_{2}$, then, when marginalizing over the randomness, we have
a new objective

${\rm arg\;min}_{W}\;\mathbb{E}_{D\sim{\rm Ber}(p)}\|Y-W(D\star
X)\|^{2}_{2}\,,$

This is equivalent to, with $\Gamma=({\rm diag}(X^{\top}X))^{\frac{1}{2}}$,

${\rm arg\;min}_{W}\;\|Y-pWX\|^{2}_{2}+p(1-p)\|\Gamma W\|^{2}_{2}\,,$

Hence, this is equivalent to a Bayesian ridge regression with a $g$-prior as
an objective function and reduces the likelihood of over-reliance on small
sets of input data in training.

##  3 Deep Learning via Partial Least Squares

Partial least squares algorithm finds projections of the input and output
vectors $X=TP+F$ and $Y=UQ+E$ in such a way that correlation between the
projected input and output is maximized. Our DL-PLS model will introduce
nonlinearity $U=G(T)$ by assuming that $U$ is a deep learner of $T$. From
Brillinger’s result and Theorem 1 we see that linear PLS calculates $T$ and
$P$ for arbitrary $G_{L}$.

$\displaystyle Y$ $\displaystyle=UQ+E$ $\displaystyle U$ $\displaystyle=G(T)$
$\displaystyle T$ $\displaystyle=XP^{T}$

where $G$ is a deep learner. Here $X=TP+F$ is inverted to $T=XP^{T}$ as
$P^{T}P=I$.

Figure 1: Relationship between output-input $Y$ and $X$
$Y=UQ+E,\;\;U=G(T),T=XP^{T}$

Figure 1 shows the relationship between $X$ and $Y$, via the link between
score vectors. $t_{1}$ and $u_{1}$ are linear combinations of $X$ and $Y$
respectively (left and right subplots). The relationship between $t_{1}$ and
$u_{1}$ is nonlinear (middle subplot).

Although we use a composite model DL-PLS, our estimation procedure is two-
step. We first estimate the score matrices and then estimate parameters of the
deep learning function $G$. This two step process is motivated by the
Brillinger (2012). The results of Brillinger guarantee that matrices $P,Q$ are
invariant (up to proportionality) to nonlinearity, even when the true
relationship between $Y$-scores and $X$-scores is nonlinear.

There are many choices for the network architecture of the deep learner $G$.

#### PLS-ReLU

For example, $G$ is a simple feed-forward ReLU neural network, for which
sparse Bayesian priors are useful to improve the generalization (Polson and
Ročková (2018)). $U$ and $T$ are $n\times L$ matrices,

$\displaystyle T$ $\displaystyle=Z_{0}$ $\displaystyle Z_{1}$
$\displaystyle=\max(Z_{0}W_{1}+b_{1},0)$ $\displaystyle U$
$\displaystyle=Z_{1}W_{2}+b_{2}$

The weights $W_{1}$ and $W_{2}$ are to be learned.

#### PLS-Autoencoder

As the input $T$ and output $U$ have the same dimensions, one can consider
using an autoencoder network. We build on the architecture of Zhang et al.
(2019).

$\displaystyle T$ $\displaystyle=Z_{0}$ $\displaystyle Z_{1}$
$\displaystyle=g_{1}(Z_{0}W_{1}+b_{1})$ $\displaystyle U$
$\displaystyle=g_{2}(Z_{1}W_{2}+b_{2})$

where $Z_{1}$ is a $n\times l_{1}$ matrix which acts as a lower dimensional
intermediate hidden layer and $l_{1}\ll l_{0}=L$.

#### PLS-Trees

One approach is Eriksson, Trygg and Wold. The first score vector, $t_{1}=Xw$
where $w=X^{T}y/\|X^{T}y\|^{2}$. After sorting all observations from maximum
to minimum score value along $t_{1}$, they divides $X$ and $Y$ in two parts to
minimize

$(1-b)\cdot\left(a\cdot\frac{Var(y_{1})+Var(y_{2})}{Var(y)}+(1-a)\cdot\frac{Var(t_{1})+Var(t_{2})}{Var(t)}\right)+b\cdot\frac{(N_{1}-N_{2})^{2}}{(N_{1}+N_{2})^{2}}$

where $Var$’s are appropriate variances and $N$’s are sample sizes. The last
term $\frac{(N_{1}-N_{2})^{2}}{(N_{1}+N_{2})^{2}}$ acts as a regularization
penalty. The hyperparameters $a$ and $b$ are chosen by cross-validation.

#### PLS-GP

Given a new predictor matrix $X_{*}$ of size $N_{*}\times p$, the same
projection $P$ produces the corresponding $T_{*}=[t_{1,*},...,t_{L,*}]$. We
can use Gaussian process regression to predict $U_{*}=[u_{1,*},...,u_{L,*}]$
from $T_{*}$ as follows

$\displaystyle T_{*}$ $\displaystyle=X_{*}P^{T}$ $\displaystyle\hat{u}_{k,*}$
$\displaystyle=g_{k,*}(t_{k,*}),k=1,2,...,L.$

where $g_{k,*}$’s are the Gaussian process regression predictors.

$\begin{bmatrix}u\\\ u_{*}\end{bmatrix}\sim N\left(\begin{bmatrix}g\\\
g_{*}\end{bmatrix},\begin{bmatrix}K&K_{*}\\\
K_{*}^{T}&K_{**}\end{bmatrix}\right)$

where $K=K(t,t)$ is $N\times N$, $K_{*}=K(t,t_{*})$ is $N\times N_{*}$, and
$K_{**}=K(t_{*},t_{*})$ is $N_{*}\times N_{*}$. $K(\cdot,\cdot)$ is a kernel
function. The conditional mean, $g_{*}$, is given by

$\displaystyle g_{*}(t_{*})$ $\displaystyle=g(t_{*})+K_{*}^{T}K^{-1}(u-g(t)).$

Then prediction of $Y$ is

$\hat{Y}_{*}=\hat{U}_{*}Q$

Fadikar et al. (2018) use PCA to reduce the 57-dimensional output vector
together with Gaussian process regression. Alternatively, one could perform
partial least squares to reduce the dimensionality to provide a better
calibration of the simulation-based model.

To summarize,

#### Algorithm

Our overall procedure includes 5 steps

  1. 1.

Given input $X$, expand by including modeler-defined features to get an
augmented set of inputs $\phi(X)$, see Hoadley (2001).

  2. 2.

Use SVD to find eigenvalues-eigenvectors of $V=ave(\phi(X)\phi(X)^{T})$

  3. 3.

Use PLS to estimate weight matrices $Q,P$ and $Y$-scores $U$. Select
dimensionality $L$ via cross-validation and scree-plot. Let $T$ denote
$X$-scores

  4. 4.

Model $U=G(T)$ where $G$ is a deep learner. Use diagnostics plot to interpret
the model.

  5. 5.

Run MCMC on model $Y=UQ+\epsilon$ where $U=G(T)$ and $T=XP^{T}$ to obtain
predictive distribution $p(Y\mid X)$

##  4 Applications

We use two simulated and two real-world data sets to study the performance of
PLS-DL in estimating the input-output relations. In our first two simulated
examples we can see that estimated $\beta_{\mathrm{PLS}}$ provides a mechanism
to select the number of latent dimensions and to plot the input-output
relations in lower dimensional space.

###  4.1 Simulated Data: ReLU and tanh

To illustrate Brillinger’s method for identifying the activation function,
$g$, we simulate data from two neural networks with ReLU and tanh activation
functions. Asymptotically $\hat{\beta}_{OLS}$ estimates $\beta$ up to
proportionality and estimates of $g(u)$ can be obtained by plotting
$(\hat{\beta}x_{j},y_{j}),j=1,\ldots n$ and smoothing $y_{j}$ values with
$\hat{\beta}x_{j}$ near $u$.

We set $n=50$ and $\beta\in\mathbb{R}^{20}$. Figure 2 plots the OLS
$\hat{\beta}$’s versus the true $\beta$’s. We see clearly that we can estimate
the $\beta$’s up to proportionality. Figure 3 plots
$(\hat{\beta}x_{j},y_{j}),j=1,\ldots n$ and show how we can recover ReLU and
tanh fucntions albeit with noise. The scaling factor between OLS
$\hat{\beta}$’s and the true $\beta$’s is estimated together with the
activation function. That is, if $\hat{\beta}_{OLS}\approx k\beta$, then
$\hat{g}(u)\approx g(\frac{1}{k}\cdot u)$.

Figure 2: OLS $\hat{\beta}$ vs. true $\beta$: the regression model
$\hat{\beta}\sim\beta$ is shown in the top-left corner. Figure 3: $y$ vs.
$\hat{\beta}x$: the solid lines are estimated activation $\hat{g}$’s and the
dashed lines are true $g$’s.

###  4.2 Simulated Data with Independent Predictors

We start with the case when predictors are independent and PLS estimates are
the same as OLS. We start with a non-linear function

$y=\log|1+Bx|+\epsilon$

where $B=(1,2,0,0)$ and $x\sim N(0,I_{4\times 4})$ and $Var(\epsilon)=0.005$.

Figure 4: The training data set generated using the function
$y=\log|1+x^{T}B|+\epsilon$ (green) and curve fitted between first score and
$y$ using neural network (red)

We further extend this example by adding another output

$\displaystyle y_{1}=$ $\displaystyle B_{1}x+\epsilon_{2}$ $\displaystyle
y_{2}=$ $\displaystyle\log|1+B_{2}x|+\epsilon_{2}$

where $B_{1}=(0,0,2,2)$, $B_{2}=(1,2,0,0)$ , and
$Var(\epsilon)=\left(\begin{array}[]{cc}0.001&0.005\\\
0.005&0.001\end{array}\right)$

The estimated coefficients by the PLS algorithms are shown in Table 1 and we
can see that PLS estimated correctly $B_{1}$ since there is a linear relation
between $B_{1}x$ and $y$ and the estimates of $B_{2}$ are proportional to the
true values.

Table 1: Coefficients estimated by the PLS algorithm $\hat{\beta}_{1}$ | $\hat{\beta}_{2}$
---|---
$0.003$ | $0.138$
$0.001$ | $0.216$
$2.002$ | $0.034$
$2.004$ | $0.003$

Then we fitted a neural network that takes two inputs, which are first two
scores calculated by PLS and uses original $y$ as an output. The neural
network correctly learns the identity relation between the first score and
$y_{1}$ and also correctly learns the non-linear part of the relation as shown
in Figure 5.

|
---|---
(a) Training Data | (b) Fitted Data
Figure 5: Training versus PLS-DL fitted data

###  4.3 Orange Juice Data

Tenenhaus et al. (2005) applies PLS to the orange juice data. There are $q$
hedonic judgements and $p$ product characteristics. In the orange juice
dataset, we have $q=96$ judges and $p=16$ predictors. In the so-called
external analysis, a principal component analysis is carried out on $X$. Then
the $Y$ are regressed on the selected principal components after predictive
cross-validation has been used to select the number of principal components.
For example, with $L=2$. The predictive rule underlying PLS is achieved as a
superposition (as in deep learning) of the PLS components $(t_{1},t_{2})$ of
the PLS regression of $Y$ on $X$ and the coefficient vector $(p_{1},p_{2})$ of
the regression of $X$ on $(p_{1},p_{2})$. Characteristics are summarised by
$(t_{1},t_{2})$ and $\hat{Y}_{k}=c_{1k}t_{1}+c_{2k}t_{2}$.

Figure 6 shows some diagnostic plots for our DL-PLS model. The left panel
represents the reduction in unexplained variance by increasing the number of
components. This screeplot suggests that first three components can explain
almost all variation. The middle panel is a biplot for $X$ with respect to the
first two components, using SVD of $X$. Samples and variables are displayed
together, where the distances measure the similarity between
samples/variables. For example, one observes that glucose and fructose are
highly correlated. And another group of variables, including sucrose, ph1,
ph2, citric, smell.int, taste.int, odor.typi, acidity, bitter and sweet, are
also close to each other. The right panel of Figure 6 is the result of
applying Brillinger on $T$ and $U$. The nonlinear function $g$ is fitted with
a smoothing spline.

Figure 6: Left: screeplot using singular values of $X^{T}Y$. Middle: Biplot
for $X$. Right: Brillinger to find the relation between $T$ ($X$-scores) and
$U$ ($Y$-scores). $\beta$ is a diagonally dominant matrix.

Figure 7 displays $X$ and $Y$ in a correlation circle (the outer circle has
radius 1 and the inner one has radius 0.5). For each point, the coordinates in
the plot are its correlations with the first two PLS components. This plot
also indicates the correlation between variables. But unlike the biplot of $X$
in Figure 6, the variable locations are now determined by both $X$ and $Y$. We
see that the large group of variables in the $X$-biplot are now separated into
multiple subgroups. For example, sucrose, ph1, ph2, sweet, odor.typi are still
close to each other, while bitter, acidity and citric form another subgroup.

Figure 7: The correlation between $X$ ($Y$) and the first two PLS components
$t_{1},t_{2}$.

Figure 8 exhibits the scaling factor $f_{j}$ versus the number of components,
where solid lines denote PLS and dotted ones are PCR. The factors for the
first four judges are shown in the plot. For PLS, the number of components is
2. Hence the second eigen-direction is expanded and $f_{j}^{PLS}>1$
accordingly. The number of components for PCR is chosen so that the overall
shrinkage are as close as possible.

Figure 8: Scaling factor $f_{j}$ of PLS (solid lines) and PCR (dotted lines)
for the first four judges.

To predict $Y$ from its score vectors $U$ with a linear model, we assume
independent and identical normal prior, $N(0,\sigma^{2}I_{L})$ where
$\sigma^{2}=0.1$, for coefficient vectors. Using the first five observations,
we derive the corresponding posteriors as well as the posterior predictive
distributions for the last observation. Results for the first three judges are
shown in Figure 9. The parameters in the posterior predictive distribution are
shown in the upper-right. Red lines denote the actual observed values.

Figure 9: Posterior predictive distribution for the last observation of $Y$

###  4.4 Wine Data: Empirical distribution of the variables

The Wine Quality Data Set (P. Cortez, A. Cerdeira, F. Almeida, T. Matos and J.
Reis, ‘Wine Quality Data Set’, UCI Machine Learning Repository.) contains 4
898 observations with 11 features. The output wine rating is an integer
variable ranging from 0 to 10 (the observed range in the data is from 3 to 9).
The frequency of each rating is as follows:

Table 2: Frequencies of Wine Ratings rating | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---
frequency | 20 | 163 | 1457 | 2198 | 880 | 175 | 5

Most of the wine falls into ratings 5 and 6.

Figure 10 shows the histograms of 11 characteristics after log transform or
square root transform where appropriate in order for approximate Gaussianity
to hold so that we can apply Brillinger’s result. Figure 11 shows the
scatterplots of wine quality versus wine characteristics.

First we use the trick of Hoadley (2001) and expand the original 11
characteristic variables, denoted as $x_{1},x_{2},...,x_{11}$, to a much
higher dimensional set of 77 variables. These include 11 corresponding
quadratics, i.e. $x_{1}^{2},x_{2}^{2},...,x_{11}^{2}$, together with $11\times
10/2=55$ interactions, i.e. $x_{1}x_{2},...,x_{10}x_{11}$. One of the
advantages of PLS is that it can be applied in high-dimensional output and
input variables that can be highly co-linear. Hence hand-coding all possible
nonlinear functionals of characteristics and including them in the expanded
set can be performed before using SVD to find the eigen-decomposition.

Figure 10: Histograms of wine variables after transformation Figure 11: Wine
variables versus wine quality

In Table 3, we compare the in-sample adjusted $R^{2}$ and out-of-sample
predictive accuracy of three linear models, OLS, PCR and PLS. The results for
both input sets (the original 11 variables and the expanded 77 variables) are
shown. For PCR, the number of components used are 5 and 15 for the original
set and expanded set, respectively. For PLS, the number of components used are
5 and 10, respectively. As expected, the Hoadley trick improves the goodness
of fit for all three models. It also helps with the out-of-sample
classification for PCR and PLS. PLS has the best performance.

Based on the expanded variable set, we also build a neural network for quality
classification (6 classes in total) and apply Brillinger’s method. The
architecture of our feed-forward neural network is 77-16-32-16-6. The second
last layer with 16 hidden nodes is extracted as the new “features”
($\tilde{x}$). We then run 6 regressions $y_{i}$ on $\tilde{x}$ to get
$\hat{\beta}_{i}$’s, $i=1,2,...,6$. The nonlinear activation $g$ is estimated
by smooth splines. The out-of-sample classification accuracy of the neural
network is 0.889; after Brillinger’s method being used, it’s 0.885.

| Reg. with 11 variables | Reg. with 77 variables
---|---|---
| adj. $R^{2}$ | OOS acc. | adj. $R^{2}$ | OOS acc.
OLS | 0.347 | 0.59 | 0.403 | 0.575
PCR | 0.324 | 0.58 | 0.372 | 0.595
PLS | 0.347 | 0.595 | 0.389 | 0.6
Table 3: Comparison results for linear models

##  5 Discussion

By merging deep learning with partial least squares, we provide a simple
methodology for non-parametric input-out pattern matching. Partial least
squares is a powerful statistical tool that performs data reduction and
feature selection in high dimensional regression prediction. We introduce PLS-
ReLU, PLS-trees and PLS-GP. Coupled with Brillinger’s estimation insight from
Stein’s lemma and Gaussian regressors, it provides a data pre-processing tool
with wide application in statistics and machine learning.

Our procedure extends PLS by modeling the $Y$-scores ($U$) as a deep learner
($G_{L}$) of the input scores ($T$). This framework also provides a number of
diagnostic plots to aid in building deep learning predictors. For example, the
scree plot can be used to identify the number of factors. We illustrate
Brillinger’s result by showing how to recover the ReLU and tanh activation
functions in a simulation study and we analyze two classic datasets, the
orange juice preference hedonic regression (Tenenhaus et al. (2005)) and the
wine quality dataset where we use the variable expansion approach of Breiman
(2001) and Hoadley (2001).

Partial least squares can be used in many well-known machine learning models,
such as kernel regression (Rosipal and Trejo (2001), Rosipal et al. (2001)),
Gaussian processes (Gramacy and Lian (2012)), tree models. Higdon et al.
(2008) uses a PCA decomposition of the $Y$-variable before applying a
nonlinear regression method. This can be addressed using our methodology.
There are many directions for future research. For example, PLS for image
processing, and extending Theorem 1 for larger classes of deep learners.

## References

  * Antoniadis et al. (2004) Antoniadis, A., G. Grégoire, and I. W. McKeague (2004).  Bayesian estimation in single-index models.  Statistica Sinica, 1147–1164.
  * Breiman (2001) Breiman, L. (2001).  Statistical modeling: The two cultures (with comments and a rejoinder by the author).  Statistical science 16(3), 199–231.
  * Brillinger (2012) Brillinger, D. R. (2012).  A Generalized Linear Model With “Gaussian” Regressor Variables.  In P. Guttorp and D. Brillinger (Eds.), Selected Works of David Brillinger, Selected Works in Probability and Statistics, pp. 589–606. New York, NY.
  * Cook (2007) Cook, R. D. (2007).  Fisher lecture: Dimension reduction in regression.  Statistical Science 22(1), 1–26.
  * Cook (2009) Cook, R. D. (2009).  Regression graphics: Ideas for studying regressions through graphics, Volume 482.  John Wiley & Sons.
  * Cox (1968) Cox, D. R. (1968).  Notes on some aspects of regression analysis.  Journal of the Royal Statistical Society: Series A (General) 131(3), 265–279.
  * Cox and Reid (2000) Cox, D. R. and N. Reid (2000).  The theory of the design of experiments.  CRC Press.
  * Diaconis and Shahshahani (1984) Diaconis, P. and M. Shahshahani (1984).  On nonlinear functions of linear combinations.  SIAM Journal on Scientific and Statistical Computing 5(1), 175–191.
  * Erdogdu et al. (2016) Erdogdu, M. A., L. H. Dicker, and M. Bayati (2016).  Scaled least squares estimator for glms in large-scale problems.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 29, pp. 3324–3332.
  * Eriksson et al. (2009) Eriksson, L., J. Trygg, and S. Wold (2009).  Pls-trees®, a top-down clustering approach.  Journal of Chemometrics 23(11), 569–580.
  * Fadikar et al. (2018) Fadikar, A., D. Higdon, J. Chen, B. Lewis, S. Venkatramanan, and M. Marathe (2018, Jan).  Calibrating a stochastic, agent-based model using quantile-based emulation.  SIAM/ASA Journal on Uncertainty Quantification 6(4), 1685–1706.
  * Fearn (1983) Fearn, T. (1983).  A Misuse of Ridge Regression in the Calibration of a Near Infrared Reflectance Instrument.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 32(1), 73–79.
  * Fisher (1922) Fisher, R. A. (1922).  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 222(594-604), 309–368.
  * Frank (1994) Frank, I. (1994).  NNPPSS: Neural networks based on pcr and pls components nonlinearized by smoothers and splines.  Technical Report TDL 92-03, BIRL.
  * Frank and Friedman (1993) Frank, I. E. and J. H. Friedman (1993).  A Statistical View of Some Chemometrics Regression Tools.  Technometrics 35(2), 109–135.
  * Friedman and Stuetzle (1981) Friedman, J. H. and W. Stuetzle (1981).  Projection pursuit regression.  Journal of the American statistical Association 76(376), 817–823.
  * Golub and Van Loan (2013) Golub, G. H. and C. F. Van Loan (2013).  Matrix computations.  JHU press.
  * Gramacy and Lee (2008) Gramacy, R. B. and H. K. H. Lee (2008).  Bayesian treed gaussian process models with an application to computer modeling.  Journal of the American Statistical Association 103(483), 1119–1130.
  * Gramacy and Lian (2012) Gramacy, R. B. and H. Lian (2012).  Gaussian process single-index models as emulators for computer experiments.  Technometrics 54(1), 30–41.
  * Gron et al. (2012) Gron, A., B. N. Jørgensen, and N. G. Polson (2012).  Optimal portfolio choice and stochastic volatility.  Applied Stochastic Models in Business and Industry 28(1), 1–15.
  * Hahn et al. (2013) Hahn, P. R., C. M. Carvalho, and S. Mukherjee (2013).  Partial Factor Modeling: Predictor-Dependent Shrinkage for Linear Regression.  Journal of the American Statistical Association 108(503), 999–1008.
  * Hall and Li (1993) Hall, P. and K.-C. Li (1993, 06).  On almost linearity of low dimensional projections from high dimensional data.  Ann. Statist. 21(2), 867–889.
  * Haykin (1998) Haykin, S. (1998).  Neural Networks: A Comprehensive Foundation (Subsequent edition ed.).  Upper Saddle River, N.J: Prentice Hall.
  * Helland (1988) Helland, I. S. (1988, January).  On the structure of partial least squares regression.  Communications in Statistics - Simulation and Computation 17(2), 581–607.
  * Helland (1990) Helland, I. S. (1990).  Partial Least Squares Regression and Statistical Models.  Scandinavian Journal of Statistics 17(2), 97–114.
  * Higdon et al. (2008) Higdon, D., J. Gattiker, B. Williams, and M. Rightley (2008).  Computer model calibration using high-dimensional output.  Journal of the American Statistical Association 103(482), 570–583.
  * Hoadley (2001) Hoadley, B. (2001).  Statistical Modeling: The Two Cultures: Comment.  Statistical Science 16(3), 220–224.
  * Huber (1985) Huber, P. J. (1985).  Projection pursuit.  The annals of Statistics, 435–475.
  * Jiang and Liu (2014) Jiang, B. and J. S. Liu (2014).  Variable selection for general index models via sliced inverse regression.  The Annals of Statistics 42(5), 1751 – 1786.
  * Jones (1992) Jones, L. K. (1992, 03).  A simple lemma on greedy approximation in hilbert space and convergence rates for projection pursuit regression and neural network training.  Ann. Statist. 20(1), 608–613.
  * Klusowski and Barron (2016) Klusowski, J. M. and A. R. Barron (2016).  Risk bounds for high-dimensional ridge function combinations including neural networks.  arXiv preprint arXiv:1607.01434.
  * Kramer (1991) Kramer, M. A. (1991).  Nonlinear principal component analysis using autoassociative neural networks.  AIChE journal 37(2), 233–243.
  * Li (1991) Li, K.-C. (1991).  Sliced inverse regression for dimension reduction.  Journal of the American Statistical Association 86(414), 316–327.
  * Li and Nachtsheim (2007) Li, L. and C. J. Nachtsheim (2007).  Comment: Fisher lecture: Dimension reduction in regression.  Statistical Science 22(1), 36–39.
  * Lindley (1968) Lindley, D. V. (1968).  The choice of variables in multiple regression.  Journal of the Royal Statistical Society: Series B (Methodological) 30(1), 31–53.
  * Liu et al. (2019) Liu, H., C. Yang, B. Carlsson, S. J. Qin, and C. Yoo (2019).  Dynamic nonlinear partial least squares modeling using gaussian process regression.  Industrial & Engineering Chemistry Research 58(36), 16676–16686.
  * Mallows (1973) Mallows, C. L. (1973).  Some Comments on CP.  Technometrics 15(4), 661–675.  Publisher: [Taylor & Francis, Ltd., American Statistical Association, American Society for Quality].
  * Malthouse et al. (1997) Malthouse, E. C., A. C. Tamhane, and R. Mah (1997).  Nonlinear partial least squares.  Computers & Chemical Engineering 21(8), 875–890.
  * Manne (1987) Manne, R. (1987).  Analysis of two partial-least-squares algorithms for multivariate calibration.  Chemometrics and Intelligent Laboratory Systems 2(1), 187–197.
  * Naik and Tsai (2000) Naik, P. and C.-L. Tsai (2000).  Partial least squares estimator for single-index models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 62(4), 763–771.
  * Pinkus (2015) Pinkus, A. (2015).  Ridge functions, Volume 205.  Cambridge University Press.
  * Polson and Sokolov (2020) Polson, N. and V. Sokolov (2020).  Deep learning: Computational aspects.  Wiley Interdisciplinary Reviews: Computational Statistics 12(5), e1500.
  * Polson and Ročková (2018) Polson, N. G. and V. Ročková (2018).  Posterior concentration for sparse deep learning.  Advances in Neural Information Processing Systems 31.
  * Polson and Scott (2010) Polson, N. G. and J. G. Scott (2010).  Shrink Globally, Act Locally: Sparse Bayesian Regularization and Prediction, Volume 105.  Oxford University Press.  Bayesian Statistics 9.
  * Polson and Scott (2012) Polson, N. G. and J. G. Scott (2012).  Local shrinkage rules, Lévy processes and regularized regression.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74(2), 287–311.
  * Polson and Sokolov (2017) Polson, N. G. and V. Sokolov (2017).  Deep Learning: A Bayesian Perspective.  Bayesian Analysis 12(4), 1275–1304.
  * Rosipal et al. (2001) Rosipal, R., M. Girolami, L. J. Trejo, and A. Cichocki (2001).  Kernel PCA for Feature Extraction and De-Noising in Nonlinear Regression.  Neural Computing & Applications 10(3), 231–243.
  * Rosipal and Trejo (2001) Rosipal, R. and L. J. Trejo (2001).  Kernel partial least squares regression in reproducing kernel hilbert space.  Journal of machine learning research 2(Dec), 97–123.
  * Stone and Brooks (1990) Stone, M. and R. J. Brooks (1990).  Continuum Regression: Cross-Validated Sequentially Constructed Prediction Embracing Ordinary Least Squares, Partial Least Squares and Principal Components Regression.  Journal of the Royal Statistical Society. Series B (Methodological) 52(2), 237–269.
  * Tenenhaus et al. (2005) Tenenhaus, M., J. Pages, L. Ambroisine, and C. Guinot (2005).  PLS methodology to study relationships between hedonic judgements and product characteristics.  Food quality and preference 16(4), 315–325.
  * Tran et al. (2020) Tran, M.-N., N. Nguyen, D. Nott, and R. Kohn (2020).  Bayesian Deep Net GLM and GLMM.  Journal of Computational and Graphical Statistics 29(1), 97–113.
  * Wold (1975) Wold, H. (1975).  Soft Modelling by Latent Variables: The Non-Linear Iterative Partial Least Squares (NIPALS) Approach.  Journal of Applied Probability 12(S1), 117–142.
  * Wold et al. (1984) Wold, S., A. Ruhe, H. Wold, and W. J. Dunn, III (1984).  The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS) Approach to Generalized Inverses.  SIAM Journal on Scientific and Statistical Computing 5(3), 735–743.
  * Xia et al. (2009) Xia, Y., H. Tong, W. K. Li, and L.-X. Zhu (2009).  An adaptive estimation of dimension reduction space.  In Exploration Of A Nonlinear World: An Appreciation of Howell Tong’s Contributions to Statistics, pp. 299–346. World Scientific.
  * Zhang et al. (2019) Zhang, C., T. Peng, J. Zhou, J. Ji, and X. Wang (2019).  An improved autoencoder and partial least squares regression-based extreme learning machine model for pump turbine characteristics.  Applied Sciences 9(19), 3987.