# Dual-sPLS: a family of Dual Sparse Partial Least Squares regressions for feature selection and prediction with tunable sparsity; evaluation on simulated and near-infrared (NIR) data

Louna Alsouki, Laurent Duval, Clément Marteau, Rami El Haddad, François Wahl

► To cite this version:

Louna Alsouki, Laurent Duval, Clément Marteau, Rami El Haddad, François Wahl. Dual-sPLS: a family of Dual Sparse Partial Least Squares regressions for feature selection and prediction with tunable sparsity; evaluation on simulated and near-infrared (NIR) data. 2023. hal-03957532

HAL Id: hal-03957532

https://ifp.hal.science/hal-03957532

Preprint submitted on 26 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dual-sPLS: a family of Dual Sparse Partial Least Squares regressions for feature selection and prediction with tunable sparsity; evaluation on simulated and near-infrared (NIR) data

Louna Alsouki<sup>a,c</sup>, Laurent Duval<sup>b</sup>, Clément Marteau<sup>c</sup>, Rami El Haddad<sup>a</sup>,
François Wahl<sup>b,c</sup>

<sup>a</sup>Laboratoire de Mathématiques et Applications, U. R. Mathématiques et modélisation, Faculté des sciences, Université Saint-Joseph, B.P. 7-5208, Mar Mikhaël Beyrouth, 1104 2020, Liban

<sup>b</sup>IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, Rueil-Malmaison, 92852, France

<sup>c</sup>Institut Camille Jordan, Université Claude Bernard Lyon 1, 43 boulevard du 11 novembre 1918, Villeurbanne, 69100, France

## Abstract

Relating a set of variables **X** to a response **y** is crucial in chemometrics. A quantitative prediction objective can be enriched by qualitative data interpretation, for instance by locating the most influential features. When high-dimensional problems arise, dimension reduction techniques can be used. Most notable are projections (e.g. Partial Least Squares or PLS) or variable selections (e.g. lasso). Sparse partial least squares combine both strategies, by blending variable selection into PLS. The variant presented in this paper, Dual-sPLS, generalizes the classical PLS1 algorithm. It provides balance between accurate prediction and efficient interpretation. It is based on penalizations inspired by classical regression methods (lasso, group lasso, least squares, ridge) and uses the dual norm notion. The resulting sparsity is enforced by an intuitive shrinking ratio parameter. Dual-sPLS favorably compares to similar regression methods, on simulated and real chemical data. Code is provided as an open-source package in R: https://CRAN.R-project.org/package=dual.spls.

**Keywords:** Partial least squares, lasso, ridge, regression, sparsity, dual norm, chemometrics, machine learning
# 1. Introduction

Two main features of chemometrics reside in first, providing reliable inference and second, offering interpretability of chemical data sources. On the one hand, one may expect to estimate, within a given precision, responses $Y \in \mathbb{R}^{N \times Q}$ (e.g. hydrocarbon properties: viscosity, density, cetane number [1]) from spectra or variables represented by quantities $X \in \mathbb{R}^{N \times P}$ (nuclear magnitude resonance or NMR, chromatography, infrared spectroscopy, etc. [2]). It aims at relating a target $Y$ to $X$ through a predictive model: for instance, NMR spectra can be linked to viscosity with predictive purposes. On the other hand, one also wishes to interpret how variables in $X$ influence quantities $Y$: for instance, which spectral bands in NMR affect most viscosity index (see e.g. [3])? This can be transcribed by a regression model, often considered linear:

$$y = X\beta + \epsilon, \quad (1)$$

where $\epsilon$ is expected to be independent of $X$, with zero mean. With the growing size of consolidated analytical chemistry databases, chemometrics still require methodologies to 1) provide accurate predictions 2) extract pertinent knowledge or offer useful insights on measurements 3) combine heterogeneous or high-dimensional data sources. When the number $P$ of variables (samples) is far greater than the number of observations (signals) $N$ ($P \gg N$), naive statistical models risk overfitting. This notably happens in standard least squares optimizations. Dimension reduction techniques are generic approaches to deal with high dimensionality. They include projection methods or variable selection algorithms. Commonly used strategies start with PCA/PCR (principal component analysis/regression), performed only on explanatory variables in $X$. They however do not incorporate information held by the response $Y$. Partial least squares (PLS) [4, 5], also called projection onto latent structures, is therefore common in chemometrics, with better prediction-prone latent components. However, PLS sometimes lacks appropriate interpretability.

As for variable selection, one often resorts to the lasso algorithm (least absolute shrinkage and selection operator [6]). Shrinkage induces a form of sparsity, which amounts to selecting important variables. It is however known to be sensitive to data types. It does not always yield interpretable coefficients. Blends of the two above – dimension reduction and variable selection – have recent avatars called sparse PLS (sPLS). While they enforce lower dimensional decompositions, they do not always provide chemically perti-
nent feature localization for physico-analytical measurements. Thereby, we propose a dual sparse PLS family dedicated to one dimensional or univariate responses: $\mathbf{y} = \mathbf{Y}$, with $\mathbf{y} \in \mathbb{R}^N$. It generalizes the standard PLS1 algorithm by supplementing it with adequate penalties. This formally provides a unified formulation for regression methods in the spirit of the lasso mentioned above, and also least-squares or ridge, all blended in a PLS formalism. It also allows *variable grouping*: the possibility to gather explanatory variables into more meaningful subsets (contiguous samples around a peak, disjoint spectral bands associated to a compound). This can be used to combine different physico-chemical modalities. Resolution resorts to the dual norm of the chosen Dual-sPLS penalty. This new method has many advantages:

1. predictions match or outperform state-of-the-art or comparable methods,

2. in the different norm options we considered, they additionally yield sparse representations of both simulated and real chemical near infrared data, even singular, a frequent ill-conditioning issue in high dimension,

3. they finally offer a interpretable localization of features.

Those three properties combined offer alternative surrogates to classical approaches (PLS, lasso, least squares, ridge). It permits both accurate inference and pertinent domain-related interpretation.

The paper is structured as follows: setting notations, we briefly revise in Section 2 the background of the PLS, recall classical variable selection methods and evoke their blending in sparse PLS schemes, previously proposed. Then, in Section 3, we explain principles behind the Dual-sPLS family and detail the list of norm penalties and their algorithms in three main instances: the (group) lasso form —being the most important— and least squares and ridge forms. Thereafter Section 4 explicits tested data (simulated and real) and the choices of model settings, calibration and validation. Each of the three penalties types are extensively benchmarked in Section 5. We finally draw concluding remarks with perspectives in Section 6 and supplementary material in the appendix.

## Notation and definitions

Matrices, vectors and scalars are denoted by boldface uppercase letters, boldface lowercase and light lowercase letters respectively, e.g. $\mathbf{X}$, $\mathbf{y}$ and $\mathbf{\lambda}$.
The transpose of matrix $\mathbf{X}$ is $\mathbf{X}^T$. The identity matrix of size $P$ is represented by $I_P$. The $\ell_1$-norm and the $\ell_2$-norm of vector a $\mathbf{w}$ of length $P$ are

$$ \|\mathbf{w}\|_1 = \sum_{p=1}^{P} |w_p| \quad \text{and} \quad \|\mathbf{w}\|_2 = \sqrt{\sum_{p=1}^{P} |w_p|^2}. \qquad (2) $$

We denote by $\ell_0(\mathbf{w})$ the sparsity index or count measure [7] of the non-zero coordinates of $\mathbf{w}$ and $\ell_0^c(\mathbf{w})$ its complement i.e. $\ell_0^c(\mathbf{w}) = P - \ell_0(\mathbf{w})$. To choose the number of latent variables we rely on the mean squared error (MSE) expressed as

$$ \text{MSE} = \frac{1}{N} \sum_{n=1}^{N} (y_n - \hat{y}_n)^2, \qquad (3) $$

for a response vector $\mathbf{y}$ of $N$ observations and a given estimate $\hat{\mathbf{y}}$. For performance evaluation, we choose the root mean squares error (RMSE), the mean absolute error (MAE) and the determination coefficient ($R^2$):

$$ \text{RMSE} = \sqrt{\frac{1}{N} \sum_{n=1}^{N} (y_n - \hat{y}_n)^2} = \frac{1}{\sqrt{N}} \|\mathbf{y} - \hat{\mathbf{y}}\|_2, \qquad (4) $$

$$ \text{MAE} = \frac{1}{N} \sum_{n=1}^{N} |y_n - \hat{y}_n| = \frac{1}{N} \|\mathbf{y} - \hat{\mathbf{y}}\|_1, \qquad (5) $$

$$ R^2 = \frac{\sum_{n=1}^{N} (y_n - \bar{y})^2}{\sum_{n=1}^{N} (y_n - \hat{y}_n)^2} \quad \text{where} \quad \bar{y} = \frac{\sum_{n=1}^{N} y_n}{N}. \qquad (6) $$

The vector of signs of $\mathbf{w}$ is noted $\text{sign}(\mathbf{w})$, and $(\mathbf{w})_+$ is the vector composed of $\mathbf{w}_p$ if $\mathbf{w}_p \ge 0$ and $0$ if $\mathbf{w}_p < 0$<fnref n="1" />.

In the following, we assume that the matrix $\mathbf{X} \in \mathbb{R}^{N \times P}$ of independent variables and the response vector $\mathbf{y} \in \mathbb{R}^N$ are mean-centered. We use the convention where columns denote variables and rows observations.

<footnote n="1">¹It corresponds to the Rectified Linear Unit (ReLU), a popular activation function for neural networks.</footnote>
## 2. Background

### 2.1. Partial Least Squares (PLS)

PLS originated from econometrics [8]. It was progressively and successfully applied to other fields [9]: social and behavioral sciences, biosciences from bioinformatics [10] to neuroimaging [11], and chemometrics [4, 5]. It denotes a class of methods aimed at explaining the relationship between explanatory data and responses with the help of latent variables. They boast the management of both formative and reflective measurements, require low sample sizes and mild distributional assumptions.

PLS avatars root on projecting data onto a lower dimensional space using new orthogonal directions constructed as linear combinations of original variables. Its principle consists in compressing the predictor **X** into a smaller score matrix of *M* < *P* variables. When Principal Component Analysis (PCA) [12] ought to best summarize **X** by taking into account only the correlation between the variables in **X**, the PLS steps up and also consider the covariance between **X** and **y**. The latent space is spanned by *M* components written as **t**<sub>*m*</sub> = **X****w**<sub>*m*</sub>, *m* ∈ {1, ..., *M*}. Weights **w**<sub>*m*</sub> are constructed in order to obtain an orthogonal basis. Several algorithms have been proposed. NIPALS (non-linear iterative partial least squares) [13] and SIMPLS [14] are most popular. When applied to a one-dimensional response, as in our case, both are shown to be equivalent. They solve the following optimization problem for the first component:

$$ \max_{\mathbf{w}} (\mathbf{y}^T \mathbf{X} \mathbf{w}) \quad \text{s.t.} \quad \|\mathbf{w}\|_2 = 1. \qquad (7) $$

The convex Problem (7) can be solved with Lagrange multipliers. For $\mu > 0$, it rewrites:

$$ \min_{\mathbf{w}} L(\mathbf{w}) \quad \text{where} \quad L(\mathbf{w}) = -\mathbf{z}^T \mathbf{w} + \mu(\|\mathbf{w}\|_2 - 1) \quad \text{and} \quad \mathbf{z} = \mathbf{X}^T \mathbf{y}. \qquad (8) $$

Solving (8) leads to

$$ \mathbf{w} = \mathbf{X}^T \mathbf{y}. \qquad (9) $$

The PLS algorithm uses the weight vector **w** to project regressor **X** into score vector **t** = **X****w**. NIPALS iteratively computes weight vectors by deflation while SIMPLS is more direct. Let $\mathcal{P}_{\mathbf{t}_{m-1}}$ denotes the orthogonal projection onto the space spanned by components $\mathbf{t}_1, \dots, \mathbf{t}_{m-1}$. The algorithm considers the part of **X** that is orthogonal to $\mathbf{t}_k, k < m$. For the $m^{th}$ component, **X** is replaced by $\mathbf{X}_m$ such that:

$$ \mathbf{X}_m = \mathbf{X} - \mathcal{P}_{\mathbf{t}_1, \dots, \mathbf{t}_{m-1}} \mathbf{X} = \mathbf{X}_{m-1} - \mathcal{P}_{\mathbf{t}_{m-1}} \mathbf{X}_{m-1}. \qquad (10) $$
The PLS1 algorithm is described in Algorithm 1.

<caption>**Algorithm 1 PLS1**</caption>

Input: $\mathbf{X}, \mathbf{y}, M$
$\mathbf{X}_1 = \mathbf{X}$
**for** $m = 1, \dots, M$ **do**
&emsp; $\mathbf{w}_m = \mathbf{X}_m^T \mathbf{y}$ (weight vector computation)
&emsp; $\mathbf{t}_m = \mathbf{X}_m \mathbf{w}_m$ (component construction)
&emsp; $\mathbf{X}_{m+1} = \mathbf{X}_m - \mathcal{P}_{\mathbf{t}_m} \mathbf{X}_m$ (deflation)
**end for**

PLS thus projects $\mathbf{X}$ onto the space of lower dimension spanned by the loadings $\mathbf{w}_1, \dots, \mathbf{w}_m$, in order to replace $\mathbf{X}$ by $\mathbf{T}=\mathbf{X}\mathbf{W}$, where $\mathbf{T} \in \mathbb{R}^{N \times M}$. The PLS regression fitted values for $M$ components is given by:

$$ \hat{\mathbf{y}} = \mathbf{T}\hat{\beta} = \mathbf{T}(\mathbf{T}^T\mathbf{T})^{-1}\mathbf{T}\mathbf{y} = \mathbf{X}\mathbf{W}(\mathbf{T}^T\mathbf{T})^{-1}\mathbf{T}\mathbf{y}. \quad (11) $$

Based on the above, PLS regression coefficients are computed as:

$$ \hat{\beta}^{PLS} = \mathbf{W}(\mathbf{T}^T\mathbf{T})^{-1}\mathbf{T}\mathbf{y}. \quad (12) $$

## 2.2. Least absolute shrinkage and selection operator

By selecting the most important features, variable selection produces a less complicated model. It has the potential advantage of being easier to handle than the complete full set of variables. The optimization problem in standard linear regression is stated as:

$$ \arg\min_{\beta} \|\mathbf{y} - \mathbf{X}\beta\|_2^2. \quad (13) $$

Provided $\mathbf{X}$ has full column rank, the ordinary least squares (LS) estimation is $\hat{\mathbf{y}}_{LS} = \mathcal{P}_{[\mathbf{x}]\mathbf{y}}\mathbf{y}$, where $[\mathbf{x}]$ is the space spanned by the columns of $\mathbf{X}$. In other terms, $\hat{\beta}^{LS} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$. A popular sparsity-based approach is the lasso developed by Tibshirani in 1996 [6]. It is reknown for its $\ell_1$ penalty scheme that shrinks less relevant variables to zero. It is obtained by solving:

$$ \arg\min_{\beta} \|\mathbf{y} - \mathbf{X}\beta\|_2^2 \quad \text{subject to} \quad \|\beta\|_1 \le \lambda. \quad (14) $$

Threshold parameter $\lambda > 0$ controls the extent of shrinkage applied to the estimate; that is, the number $\ell_0^c$ of coefficients set to zero. An appropriate
$\lambda$ is important to get interpretable results. If $\hat{\beta}^{LS}$ exists, as mentioned in [6], then for a $\lambda \ge \|\hat{\beta}^{LS}\|_1$, the lasso estimate $\hat{\beta}^1$ is equal to the ordinary least square solution. And for $\lambda = \frac{\|\hat{\beta}^{LS}\|_1}{2}$, it selects on average half of the variables. We can reformulate (14) as

$$ \arg\min_{\beta} \frac{1}{2} \|\mathbf{y} - \mathbf{X}\beta\|_2^2 + t\|\beta\|_1. \qquad (15) $$

Note that there is a (non-explicit) correspondence between parameters $\lambda$ and $t$. In the orthonormal design case, i.e. $\mathbf{X}^T\mathbf{X} = I_P$, there exists $\hat{\beta}^1$ closed form solution called *soft thresholding* verifying

$$ \hat{\beta}_p^1 = \text{sign}(\hat{\beta}_p^{\text{LS}})(|\hat{\beta}_p^{\text{LS}}| - \lambda)_+ \quad \forall p \in \{1, \dots, P\}. \qquad (16) $$

Coefficients whose magnitude is smaller than $\lambda$ are set to zero. Amplitudes of the others are shrunk with respect of the threshold. While proved successful for numerous applications, some drawbacks are reported [15, 16]. Some are: 1) non strict convexity of the criterion when the number of predictors exceeds the number of observations ($P > N$) 2) algorithm saturation when $N$ variables have been selected 3) with highly correlated variables, tendency to pick mildly representative ones.

Another shrinking method is ridge regression [17] with optimization problem:

$$ \arg\min_{\beta} \frac{1}{2} \|\mathbf{y} - \mathbf{X}\beta\|_2^2 + t\|\beta\|_2. \qquad (17) $$

Its trick is to add a diagonal matrix to $(\mathbf{X}^T\mathbf{X})$ in order to overcome the singularity problem. Therefore, the solution always exists, expressed as:

$$ \hat{\beta}^r = (\mathbf{X}^T\mathbf{X} + tI_P)^{-1}\mathbf{X}^T\mathbf{y}. \qquad (18) $$

Compared to the lasso, it uses an $\ell_2$-norm instead of the $\ell_1$ penalization but retains most variables by design.

## 2.3. Blending methods: sparse Partial Least Squares (sPLS)

Sparse Partial Least Squares (sPLS) denotes a body of works adding a variable selection flavor to the standard PLS framework. We focus here on
ones specifically using lasso inspired penalties. An $l_1$-norm can be incorporated in optimization problem (7). Noting

$$ \widehat{\mathrm{Cov}}(\mathbf{X}\mathbf{w}, \mathbf{y}) = \frac{1}{N}\mathbf{w}^T\mathbf{z}, \quad \text{with} \quad \mathbf{z} = \mathbf{X}^T\mathbf{y} = N\widehat{\mathrm{Cov}}(\mathbf{X}, \mathbf{y}), \qquad (19) $$

adding the coupling parameter $\lambda_s > 0$ and orthogonality constraint on components, we get, for the first component:

$$ \hat{\mathbf{w}} = \arg \min_{\mathbf{w} \in \mathbb{R}^p} \{-\widehat{\mathrm{Cov}}(\mathbf{X}\mathbf{w}, \mathbf{y}) + \lambda_s \|\mathbf{w}\|_1\}, \quad \text{for } \mathbf{w}^T\mathbf{w} = 1. \qquad (20) $$

Problem (20) is tackled in 2008 [18] using sparse PCA [19]. Then iterative PLS [20] is combined to singular value decomposition. We denote it as sPLS$_{LeCao}$ after the first author. In 2010 [21], Problem (20) is reframed by imposing the $l_1$ penalty on a surrogate direction close to the original vector $\mathbf{w}$, providing an approximate solution with sPLS$_{Chun}$. In 2018, [22] reformulates Problem (20) using recent results from proximal optimization [23] with sPLS$_{Durif}$. In this last case, sPLS$_{Durif}$ provides an exact and closed-form solution reminiscing the soft threshold operator. Moreover, they suggest an adaptive method for computing the sPLS weight vectors using classical PLS ones.

Along the lines of methods presented above, Dual-sPLS aims at inference and interpretability: accurate predictions combined with sparse localization features for better chemometrics performance. Following [22], we also wish to provide a means to tuning the relative sparsity of the outcome. Finally, as different analytical chemistry modalities provide different insights on chemical mixtures, the Dual-sPLS is designed to naturally allow the combination of heterogeneous datasets as a byproduct of the versatile dual norm approach<fnref n="2" />.

# 3. Dual Sparse Partial Least Squares (Dual-sPLS)

## 3.1. Motivation and purposes

The classical data fidelity $l_2$-norm is often used in penalties applications. Arbitrary norm choices may not lead to trackable algorithms. However, the concept of dual norm is a means to formulate a unifying optimization framework with practical algorithmic properties.

<footnote n="2">²Application of this extension is not performed here and is subject to a later work.</footnote>
**Definition 3.1** Let $\Omega(\cdot)$ be a norm on $\mathbb{R}^p$. For any $z \in \mathbb{R}^p$, the associated dual norm, denoted $\Omega^*(\cdot)$, is defined as

$$ \Omega^*(z) = \max_{w} (z^T w) \quad \text{s.t.} \quad \Omega(w) = 1. \qquad (21) $$

Comparing (7) and (21), we find that the optimization of the PLS objective function amounts to finding the vector $w_1$ that fits the dual norm of the $l_2$-norm of $z$, where $z = X^T y$. This gives us the incentive to evaluate different norm expressions that could be used as domain-related penalizations. Thus, for any norm $\Omega(\cdot)$ used, the first component will be:

$$ \hat{w} = \arg \min_{w \in \mathbb{R}^p} \{-z^T w\}, \quad \text{s.t.} \quad \Omega(w) = 1. \qquad (22) $$

Although formulation is generic, we emphasize four types of norms that make practical sense when dealing with measurements typically available in chemometrics. We provide the corresponding R [24] package dual.spl s [25] with a complete description. It contains the following main functions, each of them being associated to specific penalty:

1. **Dual-sPLS<sub>1</sub>** (pseudo-lasso norm, d.spl s.lasso()). Similar to the sPLS Problem (20), an intuitive norm combines $l_2$ and $l_1$:

$$ \Omega(\mathbf{w}) = \lambda \|\mathbf{w}\|_1 + \|\mathbf{w}\|_2. \qquad (23) $$

2. **Dual-sPLS<sub>gl</sub>** (pseudo-group lasso norm, d.spl s. GL()). Inspired by group lasso [26], it combines groups of measurements. It applies pseudo-lasso to each group individually while constraining the total set. For G groups, $w_g$ represents the variables of the loading vector $\mathbf{w}$ that belongs to group $g$. The corresponding norm is formulated as:

$$ \Omega(\mathbf{w}) = \sum_{g=1}^{G} \alpha_g \|\mathbf{w}_g\|_2 + \lambda_g \|\mathbf{w}_g\|_1, \qquad (24) $$

where $\alpha_g \ge 0, \forall g \in \{1, \dots, G\}$ and $\sum_{g \in 1, \dots, G} \alpha_g = 1$.

3. **Dual-sPLS<sub>LS</sub>** (pseudo-least squares norm, d.spl s.LS()). It introduces $N_1$, a matrix of $p$ columns, and applies when $\mathbf{X}$ is not singular:

$$ \Omega(\mathbf{w}) = \lambda \|\mathbf{N}_1 \mathbf{w}\|_1 + \|\mathbf{X} \mathbf{w}\|_2. \qquad (25) $$

The classical least squares solution is recovered for $\lambda = 0$.
4. Dual-sPLS<sub>r</sub> (pseudo-ridge norm, d.spls.ridge()). It deals with cases where X is singular and resorts to a ridge-like penalization:

$$ \Omega(\mathbf{w}) = \lambda_1 \|\mathbf{w}\|_1 + \lambda_2 \|\mathbf{X}\mathbf{w}\|_2 + \|\mathbf{w}\|_2. \quad (26) $$

The construction of weight vectors $w_1, \dots, w_M$ differs in each of the four cases. It however follows similar steps as for the PLS. Starting with a reformulation of optimization Problem (21) and using Lagrange multipliers, we aim at iteratively minimizing the function $L(\mathbf{w}) = -\mathbf{z}^T\mathbf{w} + \mu(\Omega(\mathbf{w}) - 1)$, for $\mu > 0$. As some norms are not differentiable, we resort to the more generic notion of subgradient $\nabla\Omega(\mathbf{w})$ [23]. It identifies to the classical differential when it is defined. The subgradient of $L$ vanishes for

$$ \nabla\Omega(\mathbf{w}) = \frac{\mathbf{z}}{\mu}. \quad (27) $$

It is then sufficient to substitute the gradient — when it exists— of the considered norm of $\Omega(\mathbf{w})$ in (27).

We provide in the following a detailed analysis for the pseudo-lasso case of Dual-sPLS (see (23)) and some remarks for the other norms. In all cases we impose that $\mathbf{w}$ and $\mathbf{z}$ lie in the same orthant; it generalizes, in $n$ dimensions, the quadrant in the 2D plane or the octant in the 3D space. In other words, corresponding coordinates of $\mathbf{w}$ and $\mathbf{z}$ have the same sign.

### 3.2. Norm options (lasso, group lasso, least squares and ridge)

#### 3.2.1. Pseudo-lasso

We reconsider Equation (23). Let $\delta$ be the sign vector of $\mathbf{w}$ and $\mathbf{z}$. By differentiating $\Omega(\mathbf{w})$, we get

$$ \nabla\Omega(\mathbf{w}) = \lambda\delta + \frac{w}{\|w\|_2}, \quad (28) $$

and by substituting it in (27), we obtain

$$ \frac{\mathbf{w}}{\|\mathbf{w}\|_2} = \frac{\mathbf{z}}{\mu} - \lambda\delta. \quad (29) $$

The closed-form solution of the Dual-sPLS<sub>l</sub> optimization problem consists in zeroing coordinates whose magnitude is lower than the soft threshold $\lambda$ and
in reducing the others toward zero. Thus, for $\nu = \lambda\mu$ and $p \in \{1, \dots, P\}$, it
can be expressed as:

$$ \frac{w_p}{\|\mathbf{w}\|_2} = \frac{1}{\mu} \delta_p(|z_p| - \nu)_+ . \qquad (30) $$

A common issue is the choice of the appropriate shrinking parameter. Cross-Validation [27], evoked in Section 4.3, is popularly adopted in sparse regressions. We choose a more intuitive option. We obtain it adaptively, according to the proportion of variables that we would like to keep in the active set at each iteration. The procedure is illustrated in Figure 1. It represents the empirical cumulative distribution of sorted magnitudes of $|\mathbf{X}^T\mathbf{y}|$ from the real data $D_{\text{NIR}}$ described later in Section (4.2). Fixing a striking ratio $\zeta$ of expected zero coefficients (e.g. $\zeta = 80\%$), we select the threshold $\nu$ at iteration $m$ as depicted. As the cumulative distribution is non-decreasing, we choose the first $x$-axis value corresponding to ordinate 0.8.

<figure><img src="image_2.png" /><figcaption>Figure 1: Empirical cumulative distribution of the sorted magnitude of $\mathbf{z} = \mathbf{X}^T\mathbf{y}$ from real data $D_{\text{NIR}}$ to control sparsity.</figcaption></figure>

To guarantee the unit norm property for $\mathbf{w}$, we set $\mu = \|\mathbf{z}_\nu\|_2$ where $\mathbf{z}_\nu$ is the
vector of coordinates $\delta_p(|z_p| - \nu)_+$ for $p \in \{1, \dots, P\}$. Consequently,

$$ \mathbf{w} = \frac{\mu}{\nu \|\mathbf{z}_\nu\|_1 + \|\mathbf{z}_\nu\|_2^2} \mathbf{z}_\nu. $$

The rationale behind constraining the direction $\mathbf{w}$ instead of the regression coefficients $\hat{\beta}$ is their collinearity. Indeed, the estimator writes $\hat{\beta} = \mathbf{W}(\mathbf{T}^T\mathbf{T})^{-1}\mathbf{T}^T\mathbf{y}$. Being collinear, soft-thresholding $\mathbf{w}$ performs a variable selection at the same location in $\hat{\beta}$ coordinates. The pseudo-lasso Dual-sPLS is described in Algorithm 2.

**Algorithm 2 Dual-sPLS**

Input: $\mathbf{X}, \mathbf{y}, M$ (number of components desired), $\zeta$ (shrinking ratio)

$\mathbf{X}_1 = \mathbf{X}$

for $m = 1, \dots, M$ do

$\mathbf{z}_m = \mathbf{X}_m^T\mathbf{y}$ (weight vector)

Find $\nu$ adaptively according to $\zeta$

$\mathbf{z}_\nu = (\delta_p(|z_p| - \nu)_+)_p$ (applying the threshold), $p \in \{1, \dots, P\}$

$\mu = \|\mathbf{z}_\nu\|_2$ and $\lambda = \frac{\nu}{\mu}$

$\mathbf{w}_p = \frac{\|\mathbf{z}_\nu\|_2}{\nu \|\mathbf{z}_\nu\|_1 + \|\mathbf{z}_\nu\|_2^2} \mathbf{z}_\nu$ (loadings)

$\mathbf{t}_m = \mathbf{X}_m \mathbf{w}_m$ (component)

$\mathbf{X}_{m+1} = \mathbf{X}_m - \mathcal{P}_{\mathbf{t}_m} \mathbf{X}_m$ (deflation)

end for

$\hat{\beta} = \mathbf{W}(\mathbf{T}^T\mathbf{T})^{-1}\mathbf{T}^T\mathbf{y}$

Note that as long as $\mathbf{w}$ and $\mathbf{z}_\nu$ are collinear, the sparsity of the results only requires the computation of $\mathbf{w}$, up to a non-zero factor.

### 3.2.2. Pseudo-group lasso

Response $\mathbf{y}$ may be explained separately by explanatory variables of different nature with prediction models. Combining them appropriately is potentially beneficial both in predictive and interpretative powers. The same reasoning could be used to partition the dataset into groups.

Physico-chemical motivation resides in segmenting a spectrum into homogeneous bands or combining complementary modalities (e.g. IR and NMR) to predict the same property (e.g. viscosity, density). We consider $G$ groups,
and $\mathbf{z}_g$ sub-vector of $\mathbf{z}$ denotes variables belonging to group $g$. The group lasso inspired norm is expressed as in Equation (24). The closed-form solution is collinear to the vector $\mathbf{z}_{\nu_g}$. It is given by

$$ \mathbf{z}_{\nu_g} = \delta_g(|\mathbf{z}_g| - \nu_g)_+ \quad \text{and} \quad \mathbf{z}_{\nu} = (\mathbf{z}_{\nu_g})_{g \in \{1, \dots, G\}}, \qquad (31) $$

$\delta_g$ being the vector of signs of $\mathbf{w}_g$ and $\nu_g = \lambda_g \mu$ for $g \in \{1, \dots, G\}$. Each group is driven by its own threshold $\nu_g$. The latter can be obtained similarly as in Section 3.2.1. Note that this Dual-sPLS version reduces to the pseudo-lasso case when $G = 1$.

### 3.2.3. Pseudo-least squares and pseudo-ridge

The above can be generalized in many ways, by defining more versatile norm shapes, including notably weighted norms. One such possibility is $\forall \mathbf{w} \in \mathbb{R}^P$

$$ \Omega(\mathbf{w}) = \lambda_1 \|\mathbf{N}_1 \mathbf{w}\|_1 + \|\mathbf{N}_2 \mathbf{w}\|_2 + \lambda_2 \|\mathbf{w}\|_2. \qquad (32) $$

It is not easily solvable in general. However, an appropriate choice of matrices $\mathbf{N}_1$ and $\mathbf{N}_2$, and factors $\lambda_1$ and $\lambda_2$ allow us to recover the lasso and group lasso norms, but also several other already known concepts, like fused lasso, least squares or ridge. We focus here on two main situations whose optimization problem resolution can be obtained analytically. An obvious option heavily inspired by least squares regression sets $\mathbf{N}_2 = \mathbf{X}$ and $\lambda_2 = 0$. Its resolution supplements the traditional least squares problem with a more selective shrinkage akin to that of our pseudo-lasso. Namely we first note

$$ \Omega(\mathbf{w}) = \lambda \|\mathbf{N}_1 \mathbf{w}\|_1 + \|\mathbf{X} \mathbf{w}\|_2. \qquad (33) $$

Then for $\nu = \mu\lambda$ and $\delta$ the vector of signs of $\mathbf{N}_1 \mathbf{w}$ and $\mathbf{N}_1 \mathbf{z}$,

$$ \frac{\mathbf{w}}{\|\mathbf{X}\mathbf{w}\|_2} = (\mathbf{X}^T \mathbf{X})^{-1} \frac{\mathbf{z}}{\mu} - \lambda (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{N}_1^T \delta, \qquad (34) $$

where we have implicitly assumed that $\mathbf{X}$ has full rank. Consequently, we penalize $|\hat{\beta}^{\text{LS}}|$ instead of $|\mathbf{z}|$. For equation (34) to take a genuine pseudo-lasso form, it is sufficient that $\mathbf{N}_1$ verifies

$$ (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{N}_1^T \delta = \text{sign} \left( (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{z} \right). \qquad (35) $$
However, as it does not play a role in loadings' computation, it does not need to be computed explicitly. Thus, the coordinates of the simplified closed-form solution is:

$$ \frac{\mathbf{w}_p}{\|\mathbf{X}\mathbf{w}\|_2} = \frac{1}{\mu}\mathrm{sign}(\hat{\beta}_p^{\mathrm{LS}})(|\hat{\beta}_p^{\mathrm{LS}}| - \nu)_+, \qquad (36) $$

When $\mathbf{X}$ is singular, the above cannot hold. Meanwhile, this case can be addressed with a regularization inspired by the ridge [17]. By choosing $\mathbf{N}_1 = I_P$, $\mathbf{N}_2 = \lambda_2\mathbf{X}$ and $\lambda_2 = 1$, equation (32) writes

$$ \Omega(\mathbf{w}) = \lambda_1 \|\mathbf{w}\|_1 + \lambda_2 \|\mathbf{X}\mathbf{w}\|_2 + \|\mathbf{w}\|_2. \qquad (37) $$

It amounts to penalize $|\mathbf{z}_{\nu_2}|$ where $\mathbf{z}_{\nu_2} = (\nu_2\mathbf{X}^T\mathbf{X}+I_P)^{-1}$ and $\nu_2 = \lambda_2 \frac{\mathbf{w}}{\|\mathbf{X}\mathbf{w}\|_2}$, instead of $|\mathbf{z}|$ like in the pseudo-lasso. For $\nu_1 = \lambda_1\mu$, the closed-form solution is formulated as:

$$ \frac{\mathbf{w}}{\|\mathbf{w}\|_2} = \frac{1}{\mu}\delta(|\mathbf{z}_{\nu_2}| - \nu_1)_+. \qquad (38) $$

where $\delta = \mathrm{sign}(\mathbf{z}_{\nu_2}\mathbf{z})$. Adding the diagonal perturbation resolves the non-invertability of $\mathbf{X}^T\mathbf{X}$.

# 4. Simulated and real data, model settings, evaluation

## 4.1. Simulated sparse data: Gaussian mixtures $D_{\mathrm{SIM}}$ and $\bar{D}_{\mathrm{SIM}}$

For an in-depth analysis of machine learning algorithms, resorting to simulated data allows an unbiased access to ground truth. We thereby propose a parametrized model. It is thought to provide similarities with common analytical chemistry data, with all sparse parameters controlled. We choose a positively weighted mixture of $K$ Gaussians peaks with preset scale $\sigma^2$ and randomly picked amplitudes $A$ and locations $\mu$. They are summed as follows and uniformly sampled:

$$ \sum_{k=1}^{K} A_k \exp(-\frac{x - \mu_k}{2\sigma^2}). \qquad (39) $$

They are affected by a stochastic Gaussian contamination. The response vector $\mathbf{y}$ is defined by an explicit linear model composed of weighted sums of $\mathbf{X}$ values. Weights can be random or fixed quantities by range of indices.
<figure><img src="image_3.png" /><figcaption>Figure 2: $D_{SIM}$ (left) and $\bar{D}_{SIM}$ (right) simulated data.</figcaption></figure>

In this work, to evaluate Dual-sPLS in both precision and information loca-
tion, we devise a sparse additive model with only $S \ll P$ positive weights and
$P-S$ null weights. Namely, only $S$ variables are responsible in the construc-
tion of response $\mathbf{y}$. This information is especially beneficial to demonstrate
the strength of variable selection in sparse methods. Since we deal with
high-dimensional situations, we simulated $D_{SIM}$: 300 mixtures of 30 Gaus-
sians represented by 1000 variables (Figure 2 (left)). Highlighted red areas
denote variables involved in the computation of response $\mathbf{y}$. The correspond-
ing matrix of $D_{SIM}$ is singular and used in the evaluation of Dual-sPLS<sub>I</sub> and
Dual-sPLS<sub>R</sub>. Since the Dual-sPLS<sub>LS</sub> is only operational with invertible ma-
trices, we also simulated non-singular data matrix $\bar{D}_{SIM}$, 200 mixtures of
100 Gaussians represented by 50 variables. The response $\mathbf{y}$ corresponding
to $\bar{D}_{SIM}$ depends only on the first five and last twelve variables as shown in
Figure 2 (right).

4.2. *Real data: near-infrared (NIR) spectroscopy* $D_{\text{NIR}}$

In chemistry, complex mixtures of molecules are analyzed with different
physico-chemical methods. Besides, determining macroscopic properties is
important to their use.

The evaluation on real data is done using NIR spectra of hydrocarbon sam-
ples. NIR is based on the principle of absorption of radiation (infrared) by
matter [28]. Infrared radiations correspond to wavenumbers directly greater
than those of the visible light spectrum. The absorption of radiation de-
pends on chemical bonds, therefore a NIR spectrum encodes information
about the composition of the sample. We focus on the density property
which is obtained by standardized methods. The IFPEN dataset $D_{\text{NIR}}$ was
partly exposed in [29, 30]. It is available at http://www.laurent-duval.eu
and subject to a forthcoming publication [31]. It is composed of 208 samples
with 1557 variables. The corresponding matrix $\mathbf{X}$ is singular. Many chem-
ical data require adequate preprocessing: normalization, baseline removal
[32], deconvolution [7]. Here we simply apply a discrete derivative obtained
with a Savitzky-Golay smoothing filter [33] of degree 2 and length 15. It
serves as both a crude baseline filter and diversity enhancement operator
[34]. The NIR preprocessed dataset $D_{\text{NIR}}$ is represented in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: $D_{\text{NIR}}$: first derivative of the NIR spectra of 208 samples.</figcaption></figure>

4.3. *Model settings: number of latent component selection*

Selecting the appropriate number *M* of latent components is crucial when building a regression model. It balances between model complexity and prediction accuracy (degrees of freedom), preventing the risk of overfitting. This issue is especially important when using PLS and its extensions in chemometrics. The commonly used procedure is based on cross validation.
First, observations are split randomly several times into a calibration and a validation set. Second, candidate models using the calibration sets are constructed with different numbers of latent components. Third, each prediction
is evaluated on the validation set with MSE. The latter are averaged for each model. Finally, the smallest one with the lowest averaged MSE reveals an adequate number of latent components.

### 4.4. Calibration and validation

The evaluation of prediction models traditionally divides the dataset into two representative sets called calibration and validation. Three main methods are used. In the first one, observations are randomly selected. The second only considers the distribution of values in the response y [35, Stratified sampling]. A third class is known as Kennard and Stone method [36]. It optimizes relative distances between observations according to variables of X. In chemometrics, one may expect the existence of a yet unknown dependence between analytical measurements and properties. Taking both X and y values for a proper calibration and validation split would be desirable. The attempts of [37] to consider X and y in a single distance with appropriate weights is not straightforward. It is difficult to adequately weight variables that do not belong to the same space. We have recently proposed a Cal-ValXy for that purpose. It consists in dividing the dataset into subgroups according to the repartition of y and applying the Kennard and Stone to each subgroup. It is summarized in Algorithm 3, and extensively described in [38].

<caption>**Algorithm 3** Calibration and validation CalValXy</caption>

Input: $X$, $X_{type}$ (index of which set belongs each observation of $X$),
**Listecal** (number of calibration points to pick from each subset)
$G$ = mean($X$) (centroid)
$C_1$ = $\max_n \|x_G - x_n\|$, $n \in \{1, ..., N\}$ (first calibration point)
$s$ = subset where $C_1$ is located
**while** **Listecal** is not empty **do**
$s \leftarrow s + 1$
Find the minmax point $C$ in subset $s$
Remove $C$ from $X$ and **Listecal**
Store $C$ in a vector of calibration index **cal**
**end while**
# 5. Comparative evaluation and discussion

We benchmark each proposed Dual-SPLS regression flavor (respectively pseudo-lasso, least squares and ridge) against its classical counterpart, and comparable sparse SPLSs, when applicable. We follow a common procedure to state the main results. First, we split observations into calibration (80%) and validation (20%). We replaced the traditional Kennard and Stone method (KS) [36] — using explanatory variables **X** only — with CalValXy (cf. Section 4.4 and [38]). The latter incorporates the response variable **y** in the splitting and proves slightly better than KS. Comparative performance is assessed in both accuracy and quality of interpretation. For the first one, common objective metrics are root mean squared error (RMSE), mean absolute error (MAE), or determination coefficient ($R^2$) (see end of Section 1). As metrics yield similar outcomes, we only compare, in the topmost figures, RMSE values for either calibration (left) or validation (right), as we increase the number *M* of latent components from one to ten. For the second one, we assess both variable selection and localization by vertically stacking regression coefficients for each compared algorithm in the bottom figure. Results are extensively discussed on simulated and real data for Dual-SPLS<sub>1</sub>, and in less details for the least squares and ridge flavors. Complementary outcomes are provided in the supplementary materials.

## 5.1. Dual-sPLS pseudo-lasso evaluation ($D_{SIM}$, $D_{NIR}$)

Dual-sPLS<sub>1</sub> is compared to standard PLS, three alternative sparse PLS (sPLS<sub>LeCao</sub> [18], sPLS<sub>Chun</sub> [21], sPLS<sub>Durif</sub> [22]) and lasso [6]. Their respective parameters are selected by cross validation (Section 4.1). Both sPLS<sub>LeCao</sub> and Dual-sPLS<sub>1</sub> explicitly specify a sparsity parameter: the (approximate) proportion of variables $\zeta$ to be discarded ($\ell_0^c/P$). We set it here to 99%. We first evaluate Dual-sPLS<sub>1</sub> on simulated data $D_{SIM}$ (Section 4.1) in Figure 4. Top-left and right plots entail that accuracy (RMSE) globally improves as the number of latent variables *M* increases for all five PLS-related methods — in both calibration and validation. The lasso performance, independent on the number of components, is represented by the sixth dotted curve. From six to ten latent variables, all curves tend to plateau, with close RMSE values. Dual-SPLS<sub>1</sub>, sPLS<sub>Chun</sub> and PLS provide the best results (lowest curves). Thus, adding more components seems unnecessary. We choose six latent variables to compare coefficient localization. On Figure 4-bottom, we stack seven panels: original spectra (1) and the coefficients for: PLS (2), Dual-sPLS<sub>1</sub> (3),
sPLS<sub>LeCao</sub> (4), sPLS<sub>Chun</sub> (5), sPLS<sub>Durif</sub> (6), lasso (7). PLS coefficients (panel 2) match the shape of the simulated data (panel 1). However, it fails to localize the most important variables, unlike sparse PLS. The $\ell_0$ criterion (Section 1) quantifies the sparsity induced by each method. Dual-sPLS<sub>1</sub>, sPLS<sub>LeCao</sub> and lasso perform best, selecting as expected a small number of variables, with an $\ell_0$ value around 40 to 60. It however is not sufficient to hint at improvements in interpretability. Looking only at variables affecting the response (shaded red background in panel 1), most compared methods exhibit significant coefficients in many (useless) areas (transparent background). Only Dual-sPLS<sub>1</sub>, sPLS<sub>LeCao</sub> present concentrated coefficients that can help chemical interpretation. On this rudimentary yet explainable model, we hint that Dual-sPLS<sub>1</sub> provides a predictive quality comparable to its challengers, and is the best in providing at the same time accurate localization on simulated data, with a verifiable (yet simplified) prediction model. We are now able to evaluate the performance of Dual-sPLS<sub>1</sub> on real near-infrared data $D_{NIR}$ (4.2) for density prediction. Similarly to $D_{SIM}$, RMSE curves in Figure 5 for calibration (top-left) and validation (top-right) globally decrease with an increasing number of components. Errors plateau after six components, indicating that additional latent structure orders might be weakly helpful. The performance gap for sPLS<sub>Durif</sub> could occur as it was mainly designed for classification. Again, we assess model interpretation in Figure 5 (bottom) for six latent vectors. By nature, location of the most influential features of spectra for a specific property is yet to be unveiled. One may expect that most of the meaningful variables are located in the active parts of the signal, e.g. spectral bands with relatively higher intensities, with some others possibly in quieter wavenumber ranges. On the top panel, NIR spectra are mainly active<fnref n="3" /> from 4000 cm<sup>-1</sup> to 4800 cm<sup>-1</sup> and 5500 cm<sup>-1</sup> to 6000 cm<sup>-1</sup>. Meaningful PLS coefficients are visible on a much wider support, provoking ambiguity on the identification of spectral bands related to density. All sPLS actually have smaller support, sPLS<sub>LeCao</sub> and Dual-sPLS<sub>1</sub> being the sparsest with $\ell_0$ respectively equal to 88 and 82. The first singularity of Dual-sPLS<sub>1</sub> is the contiguous and smoothness of its coefficients. By contrast, sPLS<sub>Chun</sub> and sPLS<sub>Durif</sub> coefficients location appear to be more scattered across the wavenumber axis, in non-contiguous small chunks and even

<footnote n="3">³We do not endeavour a chemical explanation here. It ought to be substanciated in forthcoming paper [31]</footnote>
isolated spikes. The second is the absence of response in the 5500 cm⁻¹ to 6000 cm⁻¹ bands³ in Dual-sPLS₁. We are not able to chemically explain the discrepancy of absence/presence results in this band. However, Dual-sPLS₁ does not need it to remain almost as accurate as its competitors.

## 5.2. Dual-sPLS pseudo-least squares evaluation ($\bar{D}_{SIM}$)

The Dual-sPLS<sub>LS</sub> requires data to be represented by a non-singular matrix **X**, as explained in Section 3.2.3. Since real data $D_{NIR}$ is singular, we use simulated data $\bar{D}_{SIM}$ presented in Section 4.1. As the number of variables in $\bar{D}_{SIM}$ is already small, we only shrink 60% of its variables to evaluate the Dual-sPLS<sub>LS</sub> against classical least squares. The latter is denoted by dashes, as the number of latent components is meaningless in this case.

For calibration (Figure 6 top-left) the RMSE for Dual-sPLS<sub>LS</sub> decreases mildly as the number of components increases. It approaches the least squares performance. For validation (Figure 6 top-right) Dual-sPLS<sub>LS</sub> performs similarly or better than least squares all over model orders. This contrast in performance might be explained by a tendency to overfit for least squares. A better prediction performance is expected with our model. Similarly to the Dual-sPLS₁, we also choose to evaluate it with six components in the bottom of Figure 6. Again, redish regions indicate active variables for the unknown linear model. We observe an overall similarity in the dynamics of both regression coefficients: strong amplitude in the first five and last ten variables corresponding to active regions. The main difference resides in the intermediate part, irrelevant to the response. Least squares as expected shrinks inactive variables towards zero but not as much as Dual-sPLS<sub>LS</sub> does. This is exemplified in the zoomed panels, where Dual-sPLS<sub>LS</sub> exhibit much less non-zero coefficients.

## 5.3. Dual-sPLS pseudo-ridge evaluation ($D_{SIM}$, $D_{NIR}$)

Dual-sPLS<sub>r</sub> is compared to classical ridge regression (Section 2.2) either applied to simulated data $D_{SIM}$ or real data $D_{NIR}$. Ridge hyper parameter $t$ (equation (17)) is fixed using cross validation. We set $\lambda_2$ for Dual-sPLS<sub>r</sub> (equation (26)) to $\frac{1}{t}$ for easier comparison. All other parameters are kept as for Dual-sPLS₁ (Section 5.1). Looking at top-left and -right in Figure 7 Dual-sPLS<sub>r</sub> reaches a plateau for $D_{SIM}$ after five latent components. Moreover, its RMSE values are slightly lower than ridge's for both calibration and
<figure><img src="image_5.png" /><figcaption>Figure 4: Dual-sPLS₁ evaluation on simulated data *D*<sub>SIM</sub>. (Top) RMSE values for calibration (left) and validation (right) with respect to the number of latent components. (Bottom) From top to bottom: simulated data *D*<sub>SIM</sub>, regression coefficients of PLS, Dual-sPLS₁, sPLS<sub>LeCao</sub>, sPLS<sub>Chun</sub>, sPLS<sub>Durif</sub> for six components, and lasso.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Dual-sPLS₁ evaluation on real data *D*<sub>NIR</sub>. (Top) RMSE values for calibration (left) and validation (right) with respect to the number of latent components. (Bottom) From top to bottom: original data *D*<sub>NIR</sub>, regression coefficients of PLS, Dual-sPLS₁, sPLS<sub>LeCao</sub>, sPLS<sub>Chun</sub>, sPLS<sub>Durif</sub> for six components, and lasso.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Dual-sPLS<sub>LS</sub> evaluation on simulated data D̄<sub>SIM</sub>. (Top) RMSE values for calibration (left) and validation (right) with respect to the number of latent components. (Bottom) From top to bottom: simulated data D̄<sub>SIM</sub>, regression coefficients of least squares and Dual-sPLS<sub>LS</sub> for five components.</figcaption></figure>
validation. We can safely select six latent components as before. Reference coefficients for ridge are misleading because the largest ones do not reside in influencing areas. They therefore can not be used for data interpretation. By selecting only fifty variables, located in red regions governing the model, Dual-sPLS<sub>r</sub> better succeeds in both prediction and localization. Similar conclusions can be drawn for real data *D*<sub>NIR</sub> on RMSE values. Dual-sPLS<sub>r</sub> even better predicts the response *y* with only four components. Regression coefficients (Figure 8 bottom) yield comments akin to above. While ridge apparently emphasizes unimportant features, Dual-sPLS<sub>r</sub> seems more reliable in identifying of relevant variables to predict density using chemical data.

# 6. Conclusion and perspectives

We propose a family of dual sparse Partial Least Squares algorithms that broadens the compass of standard PLS. Along with competitive prediction accuracy with respect to PLS as used in chemometrics, we expect additional benefits in dimension reduction or model interpretability. This is achieved by supplementing the traditional optimization problem with well-chosen dual norms.

We chiefly validate this approach by borrowing three classical regression penalties: lasso, least-squares, ridge. Each proposed Dual-sPLS draws close to the reference in calibration/validation performance with a reduced number of latent components. This is assessed in a benchmark on both realistic simulated models and real near infrared spectroscopy data, against a standard baseline and sparse contenders. Coefficients are sieved with a user-defined sparsity target. They are well-located in influential data ranges, suggesting a means for better interpretability of the trained prediction reduced model. Pseudo-lasso and ridge Dual-sPLS avatars exhibit close colocation of selected features in both datasets despite different penalties. This suggests a robust identification of meaningful information in signals.

The Dual-sPLS framework is thus a good candidate for a host of applications. We provide it as an open-source package in R [25]. It can be prolonged to other field-favorite penalties, for instance elastic net. We plan to evaluate the alluded “pseudo-group lasso” option, to refine feature selection on important contiguous areas, or to combine datasets providing complementary information on the predicted response. To improve prediction robustness or reduce the number of necessary latent components (toward three or four instead of six), we explore additional diversity enhancement preprocessing, such as
<figure><img src="image_8.png" /><figcaption>Figure 7: Dual-sPLS<sub>r</sub> evaluation on simulated data *D*<sub>SIM</sub>. (Top) RMSE values for calibration (left) and validation (right) with respect to the number of latent components. (Bottom) From top to bottom: original data *D*<sub>SIM</sub>, regression coefficients of ridge and Dual-sPLS<sub>r</sub> for five components.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Dual-sPLS<sub>r</sub> evaluation on real data *D*<sub>NIR</sub>. (Top) RMSE values for calibration (left) and validation (right) with respect to the number of latent components. (Bottom) From top to bottom: original data *D*<sub>NIR</sub>, regression coefficients of ridge and Dual-sPLS<sub>r</sub> for five components.</figcaption></figure>
[13] H. Wold, Path models with latent variables: The NIPALS approach, in: Quantitative Sociology. International Perspectives on Mathematical and Statistical Modeling, Elsevier, 1975, pp. 307–357. doi:10.1016/b978-0-12-103950-9.50017-4.

[14] S. de Jong, SIMPLS: An alternative approach to partial least squares regression 18 (3) (1993) 251–263. doi:10.1016/0169-7439(93)85002-x.

[15] H. Zou, T. Hastie, Regularization and variable selection via the elastic net, J. R. Stat. Soc. Ser. B Stat. Methodol. 67 (2) (2005) 301–320. doi:10.1111/j.1467-9868.2005.00503.x.

[16] T. Hastie, R. Tibshirani, M. Wainwright, Statistical Learning with Sparsity: The Lasso and Generalizations, CRC Press, 2015.

[17] A. E. Hoerl, R. W. Kennard, Ridge regression: Applications to nonorthogonal problems, Technometrics 12 (1) (1970) 69–82. doi:10.1080/00401706.1970.10488635.
URL http://dx.doi.org/10.1080/00401706.1970.10488635

[18] K.-A. Lê Cao, D. Rossouw, C. Robert-Granié, P. Besse, A sparse PLS for variable selection when integrating omics data, Stat. Appl. Genet. Mol. Biol. 7 (1) (2008) 35. doi:10.2202/1544-6115.1390.

[19] H. Shen, J. Z. Huang, Sparse principal component analysis via regularized low rank matrix approximation, J. Multivar. Anal. 99 (6) (2008) 1015–1034. doi:10.1016/j.jmva.2007.06.007.

[20] M. Tenenhaus, La régression PLS. Théorie et pratique, Éditions Technip, 1998.

[21] H. Chun, S. Keleş, Sparse partial least squares regression for simultaneous dimension reduction and variable selection, J. R. Stat. Soc. Ser. B Stat. Methodol. 72 (1) (2010) 3–25. doi:10.1111/j.1467-9868.2009.00723.x.
URL http://dx.doi.org/10.1111/j.1467-9868.2009.00723.x

[22] G. Durif, L. Modolo, J. Michaelsson, J. E. Mold, S. Lambert-Lacroix, F. Picard, High dimensional classification with combined adaptive sparse PLS and logistic regression, Bioinformatics 34 (3) (2018) 485–493. doi:10.1093/bioinformatics/btx571.
[23] F. Bach, R. Jenatton, J. Mairal, G. Obozinski, Optimization with sparsity-inducing penalties, Found. Trends Mach. Learn. 4 (1) (2012) 1–106. doi:10.1561/2200000015.

[24] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria (2021). URL https://www.R-project.org/

[25] L. Alsouki, F. Wahl, G. Durif, dual.spl s: Dual sparse partial least squares regression, CRAN, R package version 0.1.2 (Oct. 2022). URL https://CRAN.R-project.org/package=dual.spl s

[26] N. Simon, J. Friedman, T. Hastie, R. Tibshirani, A sparse-group lasso, J. Comp. Graph. Stat.) 22 (2) (2013) 231–245. doi:10.1080/10618600. 2012.681250.

[27] M. Stone, Cross-validatory choice and assessment of statistical predic-
tions, J. R. Stat. Soc. Ser. B Stat. Methodol. 36 (2) (1974) 111–133.
doi:10.1111/j.2517-6161.1974.tb00994.x.

[28] J. M. Chalmers, P. R. Griffiths (Eds.), Handbook of Vibrational Spectroscopy, Wiley, 2002.

[29] J. Laxalde, C. Ruckebusch, O. Devos, N. Caillol, F. Wahl, L. Duponchel, Characterisation of heavy oils using near-infrared spectroscopy: Optimisation of pre-processing methods and variable selection, Anal. Chim. Acta 705 (1-2) (2011) 227–234. doi:10.1016/j.aca.2011.05.048.

[30] J. Laxalde, Analyse des produits lourds du pétrole par spectroscopie infrarouge, Ph.D. thesis, Université de Lille 1 (2012).

[31] L. Duval, L. Alsouki, F. Wahl, J. Laxalde, N. Caillol, IFPEN near-infrared spectroscopy dataset for property prediction: 208 NIR hydrocarbon spectra and density response, PREPRINT (2023).

[32] X. Ning, I. W. Selesnick, L. Duval, Chromatogram baseline estimation and denoising using sparsity (BEADS) 139 (2014) 156–167. doi:10.1016/j.chemolab.2014.09.014.

[33] A. Savitzky, M. J. E. Golay, Smoothing and differentiation of data by simplified least squares procedures, Anal. Chem. 36 (8) (1964) 1627—1639.
[34] L. K. DeNoyer, J. G. Dodd, Smoothing and derivatives in spectroscopy. doi:10.1002/0470027320.s4501.

[35] C.-E. Särndal, B. Swensson, J. Wretman, Model Assisted Survey Sampling, Springer, 2003.

[36] R. W. Kennard, L. A. Stone, Computer aided design of experiments, Technometrics 11 (1) (1969) 137–148. doi:10.1080/00401706.1969.10490666.

[37] H. Tian, L. Zhang, M. Li, Y. Wang, D. Sheng, J. Liu, C. Wang, Weighted SPXY method for calibration set selection for composition analysis based on near-infrared spectroscopy, Infrared Phys. Technol. 95 (2018) 88–92. doi:10.1016/j.infrared.2018.10.030.

[38] L. Alsouki, L. Duval, R. El Haddad, C. Marteau, F. Wahl, CalValXy: well-balanced and stratified calibration/validation splitting using both predictors X and response y, PREPRINT (2022).

## Appendix A. Detailed resolution of Dual-sPLSs

### Appendix A.1. Dual-sPLS pseudo-group lasso

We recall Equation (24): the Dual-sPLS<sub>gl</sub> norm case applied to optimization Problem (21). Note that here

* g represents a group of P(g) index extracted from {1, ..., P};

* G represents the number of groups;

* $w_g$ represents the values of index g in the loading vector w.

We denote $z_g$ the variables of z belonging to group g. We impose $z_g$ and $w_g$ to be in the same orthant. Let $\delta_g$ be their vector of signs. By differentiating equation (24) we obtain

$$ \frac{\partial \Omega(\mathbf{w})}{\partial w_g} = \frac{\alpha_g w_g}{\|\mathbf{w}_g\|_2} + \alpha_g \lambda_g \delta_g. \quad (A.1) $$

Using Lagrange multipliers as in Section 3.1, we compare (27) to (A.1) and obtain for $g \in \{1, ..., G\}$:
$$ \frac{\mathbf{w}_g}{\|\mathbf{w}_g\|_2} = \frac{\mathbf{z}_g}{\alpha_g \mu} - \lambda_g \delta_g, \qquad (A.2) $$

which is simplified by

$$ \frac{\mathbf{w}_g}{\|\mathbf{w}_g\|_2} = \frac{1}{\mu \alpha_g} \mathbf{z}_{\nu_g}, \qquad (A.3) $$

where

$$ \mathbf{z}_{\nu_g} = \delta_g(|\mathbf{z}_g| - \nu_g)_+ \quad \text{for } g \in \{1, \dots, G\}. \qquad (A.4) $$

Here $\nu_g = \mu\alpha_g\lambda_g$ and controls the amount of variables that we would like to shrink to zero. By applying $\ell_2$-norm to (A.3), we conclude that for $g \in \{1, \dots, G\}$,

$$ \mu = \sum_{g=1}^{G} \|\mathbf{z}_{\nu_g}\|_2 \quad \text{and} \quad \alpha_g = \frac{\|\mathbf{z}_{\nu_g}\|_2}{\mu}. \qquad (A.5) $$

The term $\|\mathbf{w}_g\|_2$ is more involved. Thus, we simply use grid search. For each group $g$, ten possible values are chosen to be tested. The selection is done by detecting the maximum value of $\|\mathbf{w}_g\|_2$ for each group $g$, denoted $\|\mathbf{w}_g\|_2^{\max}$. The latter is computed by zeroing $\|\mathbf{w}_g\|_2$ for all groups $g' \neq g$ and is expressed as:

$$ \|\mathbf{w}_g\|_2^{\max} = \frac{\mu}{\Omega_g(\mathbf{z}_{\nu_g})}. \qquad (A.6) $$

Then, ten values of each group $g$ are selected inside the interval $[0, \|\mathbf{w}_g\|_2^{\max}]$. The grid search tests all the possible combinations and retains the one that allows the smallest error. We summarize the methodology with Algorithm 4.
**Algorithm 4** Dual-sPLS<sub>gl</sub> algorithm

Input: $\mathbf{X}^1, \dots, \mathbf{X}^G$, $\mathbf{y}$, $M$ (number of components desired), $\zeta$ (shrinking ratio), $\alpha_1, \dots, \alpha_g$.

for $m = 1, \dots, M$ do

$\mathbf{X}_m = (\mathbf{X}^1, \dots, \mathbf{X}^G)$ (combining data)

$\mathbf{z}_m = \mathbf{X}_m^T \mathbf{y}$ (weight vector)

Find $\nu$ adaptively according to $\zeta$ for each group separately

$\mathbf{z}_{\nu_g} = \delta_g(|\mathbf{z}_g| - \nu_g)_+$ for $g \in \{1, \dots, G\}$ (applying the threshold)

$\mu = \sum_{g=1}^G \|\mathbf{z}_{\nu_g}\|_2$

$\alpha_g = \frac{\|\mathbf{z}_{\nu_g}\|_2}{\mu}$ and $\lambda_g = \frac{\nu_g}{\alpha_g \mu}$ for $g \in \{1, \dots, G\}$

$\|\mathbf{w}_g\|_2^{\max} = \frac{\mu}{\Omega_g(\mathbf{z}_{\nu_g})}$ for $g \in \{1, \dots, G\}$

selection of the values of $\|\mathbf{w}_g\|_2$ for each group

$\mathbf{w}_g = \frac{\|\mathbf{w}_g\|_2}{\mu \alpha_g} \mathbf{z}_{\nu_g}$ for $g \in \{1, \dots, G\}$ (loadings)

$\mathbf{w}_g = (\mathbf{w}_g)_{g=1}^G$

$\mathbf{t}_m = \mathbf{X}_m \mathbf{w}_m$ (component)

$\mathbf{X}_{m+1} = \mathbf{X}_m - \mathcal{P}_{\mathbf{t}_m} \mathbf{X}_m$ (deflation)

end for

Compute $\hat{\beta}$.

Appendix A.2. Dual-sPLS pseudo-least squares

We recall Equation (25): the Dual-sPLS<sub>LS</sub> pseudo case applied to optimization Problem (21).

We impose $\mathbf{N}_1\mathbf{z}$ and $\mathbf{N}_1\mathbf{w}$ to be in the same orthant. Let $\delta_2$ be their vector of signs. By differentiating (25) we obtain

$$ \nabla \Omega(\mathbf{w}) = \lambda \mathbf{N}_1^T \delta_2 + \frac{\mathbf{X}^T \mathbf{X} \mathbf{w}}{\|\mathbf{X} \mathbf{w}\|_2}. \quad (A.7) $$

Using Lagrange multipliers as in Section 3.1, we compare (27) to (A.7) and obtain

$$ \frac{\mathbf{w}}{\|\mathbf{X}\mathbf{w}\|_2} = (\mathbf{X}^T \mathbf{X})^{-1} \frac{\mathbf{z}}{\mu} - \lambda (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{N}_1^T \delta_2, \quad (A.8) $$
imposing the invertibility of $X^T X$. We choose $N_1$ such as

$$ (X^T X)^{-1} N_1^T \delta_2 = \text{sign}((X^T X)^{-1} z). \qquad (A.9) $$

The resolution steps are similar to the ones from Dual-sPLS but instead of applying the threshold on $z$, we apply it on $(X^T X)^{-1} z$ which is exactly the classical Least Squares regression coefficients $\hat{\beta}^{LS}$. So, the simplified solution is

$$ \frac{\mathbf{w}}{\|\mathbf{X}\mathbf{w}\|_2} = \frac{1}{\mu}\text{sign}(\hat{\beta}_{LS_j})(|\hat{\beta}_{LS_j}| - \nu)_+, \qquad (A.10) $$

where $\nu$ is chosen adaptively.

For a simpler algorithm, $\|\mathbf{X}\mathbf{w}\|_2$ is not computed as it is not mandatory in this case. Additionally, $\mathbf{w}$ only depends on $\nu$ and $\hat{\beta}_{LS}$, which means $N_1$ does not intervene in the computation of the optimal solution. Thus, proving that $N_1$ exists is enough. (A.9) implies the following

$$ N_1^T \delta_2 = (X^T X)\text{sign}((X^T X)^{-1} z). \qquad (A.11) $$

Let $\mathbf{w}$ be an eigenvector of $N_1$, and $N_1'$ be such as

$$ N_1' = N_1 - \mathbf{w}^T \mathbf{w} \quad \text{and} \quad N_1' \mathbf{w} = 0. \qquad (A.12) $$

Therefore, using (A.11) we have

$$ N_1' \delta_2 = (X^T X)\text{sign}((X^T X)^{-1} z) - \mathbf{w}\mathbf{w}^T \delta_2 \quad \text{with} \quad N_1' \mathbf{w} = 0. \qquad (A.13) $$

With $N_1$ a square matrix of $P$ variables, (A.13) is a system of $P^2$ unknowns, $P$ equations and $P$ constraints. It can be verified by an infinite number of solutions.

The following algorithm reformulates the previous steps:
**Algorithm 5** DUAL-SPLS<sub>LS</sub> ALGORITHM

Input: **X**, **y**, *M* (number of components desired), *ζ* (shrinking ratio)

$\mathbf{X}_1 = \mathbf{X}$

for *m* = 1, ..., *M* do

$\mathbf{z}_m = \mathbf{X}_m^T \mathbf{y}$ (weight vector)

$\hat{\beta}_{LS} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{z}$

Find $\nu$ adaptively according to $\zeta$ and $\hat{\beta}_{LS}$

$\mathbf{z}_\nu = (\text{sign}(\hat{\beta}_{LS})(|\hat{\beta}_{LS}| - \nu)_+) \text{ (applying the threshold)}$

$\mathbf{w}_m = \frac{\mathbf{z}_\nu}{\mu}$ (loadings)

$\mathbf{w}_m = \frac{\mathbf{w}_m}{\|\mathbf{w}\|_2}$ (normalizing loadings)

$\mathbf{t}_m = \mathbf{X}_m \mathbf{w}_m$ (component)

$\mathbf{X}_{m+1} = \mathbf{X}_m - \mathcal{P}_{\mathbf{t}_m} \mathbf{X}_m$ (deflation)

end for

Compute $\hat{\beta}$.

Appendix A.3. Dual-sPLS pseudo-ridge

We recall Equation (26): the Dual-sPLS<sub>r</sub> pseudo case applied to optimization Problem (21). We impose **z** and **w** to be in the same orthant. Let **δ** be their vector of signs. By differentiating (26), we obtain

$$ \nabla \Omega(\mathbf{w}) = \lambda_1 \delta + \lambda_2 \frac{\mathbf{X}^T \mathbf{X} \mathbf{w}}{\|\mathbf{X} \mathbf{w}\|_2} + \frac{\mathbf{w}}{\|\mathbf{w}\|_2}. \quad (A.14) $$

Using Lagrange multipliers as in Section 3.1, we compare (27) to (A.14) and obtain

$$ \frac{\mathbf{w}}{\|\mathbf{w}\|_2} = \left( \nu_2 \mathbf{X}^T \mathbf{X} + I_P \right)^{-1} (\mathbf{z} - \nu_1 \mathbf{\delta}), \quad (A.15) $$

where $\nu_1 = \lambda_1 \mu$ and $\nu_2 = \lambda_2 \frac{\|\mathbf{w}\|_2}{\|\mathbf{X} \mathbf{w}\|_2}$.

In line with Dual-sPLS<sub>1</sub>, we note $\mathbf{z}_{\mathbf{X},\nu_2} = (\nu_2 \mathbf{X}^T \mathbf{X} + I_P)^{-1} \mathbf{z}$ and $\mathbf{\delta}_{\mathbf{X}}$ its vector of signs. We exhibit a solution imposing that **w** and $\mathbf{z}_{\mathbf{X},\nu_2}$ are in the same orthant, which leads to the following reformulation of (A.15):

$$ \frac{\mathbf{w}}{\|\mathbf{w}\|_2} = \frac{1}{\mu} \delta_{\mathbf{X}}(|\mathbf{z}_{\mathbf{X},\nu_2}| - \nu_1)_+. \quad (A.16) $$
The threshold $\nu_1$ is chosen with the adaptive procedure described in Section 3.2.1 and Figure 1. However, in this case, we compare $\nu_1$ to $|z_{\mathbf{X},\nu_2}|$. Since the latter is colinear to $\mathbf{z}$, the shrinkage is adequate. Denoting $z_\nu = \delta_{\mathbf{X}}(|z_{\mathbf{X},\nu_2}| - \nu_1)_+$, simple computations lead to

$$ \mu = \|z_\nu\|_2, \qquad (A.17) $$

and

$$ \mathbf{w} = \frac{\mu}{\nu_1 \|z_\nu\|_1 + \nu_2 \|\mathbf{X}z_\nu\|_2^2 + \mu^2}. \qquad (A.18) $$

It is summarized in Algorithm 6:

**Algorithm 6 DUAL-SPLS<sub>R</sub> ALGORITHM**

Input: $\mathbf{X}, \mathbf{y}, M$ (number of components desired), $\zeta$ (shrinking ratio), $\nu_2$

$\mathbf{X}_1 = \mathbf{X}$

for $m = 1, \dots, M$ do

$\mathbf{z}_m = \mathbf{X}_m^T \mathbf{y}$ (weight vector)

$\mathbf{z}_{\mathbf{X},\nu_2} = (\nu_2 \mathbf{X}^T \mathbf{X} + I_P)^{-1} \mathbf{z}$

Find $\nu$ adaptively according to $\zeta$ and $|z_{\mathbf{X},\nu_2}|$

$\delta_{\mathbf{X}}$ vector of signs of $z_{\mathbf{X},\nu_2}$

$z_\nu = \delta_{\mathbf{X}}(|z_{\mathbf{X},\nu_2}| - \nu_1)_+$ (applying the threshold)

$\mu = \|z_\nu\|_2$ and $\lambda = \frac{\nu}{\mu}$

$\mathbf{w}_m = \frac{\mu}{\nu_1 \|z_\nu\|_1 + \nu_2 \|\mathbf{X}z_\nu\|_2^2 + \mu^2}$ (loadings)

$\mathbf{t}_m = \mathbf{X}_m \mathbf{w}_m$ (component)

$\mathbf{X}_{m+1} = \mathbf{X}_m - \mathcal{P}_{\mathbf{t}_m} \mathbf{X}_m$ (deflation)

end for

Compute $\hat{\beta}$.

**Appendix B. Complementary plots**

As mentioned in Section 5, metrics MAE and $R^2$ were also computed. They support our findings based on RMSE, as they yield similar results.
<figure><img src="image_10.png" /><figcaption>Figure B.9: Dual-sPLS<sub>1</sub> evaluation on simulated data *D*<sub>SIM</sub>. MAE (top) and R<sup>2</sup> (bottom) values for calibration (left) and validation (right) with respect to the number of latent components derived from PLS, Dual-sPLS<sub>1</sub>, sPLS<sub>LeCao</sub>, sPLS<sub>Chun</sub>, sPLS<sub>Durif</sub> and lasso regressions.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure B.10: Dual-sPLS₁ evaluation on real data *D*<sub>NIR</sub>. MAE (top) and R<sup>2</sup> (bottom) values for calibration (left) and validation (right) with respect to the number of latent components derived from PLS, Dual-sPLS₁, sPLS<sub>LeCao</sub>, sPLS<sub>Chun</sub>, sPLS<sub>Durif</sub> and lasso regressions.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure B.11: Dual-sPLS<sub>LS</sub> evaluation on simulated data $\bar{D}_{SIM}$. MAE (top) and R<sup>2</sup> (bottom) values for calibration (left) and validation (right) with respect to the number of latent components derived from Dual-sPLS<sub>LS</sub> and least squares regressions.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure B.12: Dual-sPLS<sub>r</sub> evaluation on simulated data *D*<sub>SIM</sub>. MAE (top) and R<sup>2</sup> (bottom) values for calibration (left) and validation (right) with respect to the number of latent components derived from Dual-sPLS<sub>r</sub> and ridge regressions.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure B.13: Dual-sPLS<sub>r</sub> evaluation on simulated data *D*<sub>NIR</sub>. MAE (top) and R<sup>2</sup> (bottom) values for calibration (left) and validation (right) with respect to the number of latent components derived from Dual-sPLS<sub>r</sub> and ridge regressions.</figcaption></figure>