## Jointly Sparse Global SIMPLS Regression

## Tzu-Yu Liu

University of California, Berkeley, CA, USA. e-mail: tyliu@eecs.berkeley.edu

## Laura Trinchera

NEOMA Business School, Mont-Saint-Aignan, France. e-mail: laura.trinchera@neoma-bs.fr

## Arthur Tenenhaus

Sup´ elec, Gif-sur-Yvette, France. e-mail: Arthur.Tenenhaus@supelec.fr

## Dennis Wei

IBM T. J. Watson Research Center, NY, USA. e-mail: dwei@us.ibm.com

## Alfred O. Hero ∗

University of Michigan, Ann Arbor, MI, USA.

e-mail:

hero@eecs.umich.edu

Abstract: Partial least squares (PLS) regression combines dimensionality reduction and prediction using a latent variable model. Since partial least squares regression (PLS-R) does not require matrix inversion or diagonalization, it can be applied to problems with large numbers of variables. As predictor dimension increases, variable selection becomes essential to avoid over-fitting, to provide more accurate predictors and to yield more interpretable parameters. We propose a global variable selection approach that penalizes the total number of variables across all PLS components. Put another way, the proposed global penalty encourages the selected variables to be shared among the PLS components. We formulate PLS-R with joint sparsity as a variational optimization problem with objective function equal to a novel global SIMPLS criterion plus a mixed norm sparsity penalty on the weight matrix. The mixed norm sparsity penalty is the glyph[lscript] 1 norm of the glyph[lscript] 2 norm on the weights corresponding to the same variable used over all the PLS components. A novel augmented Lagrangian method is proposed to solve the optimization problem and soft thresholding for sparsity occurs naturally as part of the iterative solution. Experiments show that the modified PLS-R attains better or as good performance with many fewer selected predictor variables.

Keywords and phrases: PLS, variable selection, dimension reduction, augmented Lagrangian optimization.

∗ to whom correspondence should be addressed.

## Contents

| 1            | Introduction . . . . . . .                                                                                                                           | . . . . . . . . . . . . . . . 2                                                                                                                      |
|--------------|------------------------------------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------|
| 2            | Partial Least Squares Regression . . . . . . . . . . . . . . . . . . . . . 3                                                                         | Partial Least Squares Regression . . . . . . . . . . . . . . . . . . . . . 3                                                                         |
|              | 2.1                                                                                                                                                  | . . . . . . . . . . . . . . . 4                                                                                                                      |
|              | 2.2                                                                                                                                                  | . . . . . . . . . . . . . . . 5                                                                                                                      |
| 3            | Sparse Partial Least Squares Regression . . . . . . . . . . . . . . . . . 6                                                                          | Sparse Partial Least Squares Regression . . . . . . . . . . . . . . . . . 6                                                                          |
|              | 3.1                                                                                                                                                  | . . . . . . . . . . . . . . . 6                                                                                                                      |
|              | 3.2                                                                                                                                                  | Regression . . . . . . . . . . . . . 7                                                                                                               |
| 4            | Algorithmic Implementation for Jointly Sparse Global SIMPLS Re- gression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8     | Algorithmic Implementation for Jointly Sparse Global SIMPLS Re- gression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8     |
| 5            | Simulation Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 12                                                                        | Simulation Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 12                                                                        |
| 6            | Application 1: Chemometrics Study . . . . . . . . . . . . . . . . . . . 15                                                                           | Application 1: Chemometrics Study . . . . . . . . . . . . . . . . . . . 15                                                                           |
| 7            | Application 2: Sparse Prediction of Disease Symptoms from Gene Ex- pression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 | Application 2: Sparse Prediction of Disease Symptoms from Gene Ex- pression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 |
| 8            | Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18                                                                      | Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18                                                                      |
| References . | References .                                                                                                                                         | . . . . . . . . . . . . . . . 20                                                                                                                     |

## 1. Introduction

With advancing technology comes the need to extract information from increasingly high-dimensional data, whereas the number of samples is often limited. Dimension reduction techniques and models incorporating sparsity become important solution strategies. Partial least squares regression (PLS-R) combines dimensionality reduction and prediction using a latent variable model. It was first developed for regression analysis in chemometrics [1, 2], and has been successfully applied to many different areas, including sensory science and, more recently, genetics [3, 4, 5, 6, 7].

Moreover, PLS-R algorithm is designed precisely to operate with high dimensional data. Since the first proposed algorithm does not require matrix inversion nor diagonalization but deflation to find the latent components, it can be applied to problems with large numbers of variables. The latent components reduce the dimension by constructing linear combinations of the predictors, which successfully solved the collinearity problems in chemometrics [8]. However, the linear combinations are built on all the predictors. The resulting PLS-R model tends to overfit when the number of predictors increases for a fixed number of samples. Therefore, variable selection becomes essential for PLS-R in high-dimensional sample-limited problems. It not only avoids over-fitting, but also provides more accurate predictors and yields more interpretable estimates. For this reason sparse PLS-R was developed by H. Chun and S. Keles [9]. The sparse PLS-R algorithm performs variable selection and dimension reduction simultaneously using an glyph[lscript] 1 type variable selection penalty. However, the glyph[lscript] 1 penalty used in [9] penalizes each variable in each component independently and this can result in different sets of variables being selected for each PLS component leading to an excessively large number of variables.

In this work we propose a global criterion for PLS that changes the sequential optimization for a K component model in Statistically Inspired Modification of PLS (SIMPLS) [10] into a unified optimization formulation, which we refer to as global SIMPLS. This enables us to perform global variable selection, which penalizes the total number of variables across all PLS components. We formulate PLS-R with global sparsity as a variational optimization problem with the objective function equal to the global SIMPLS criterion plus a mixed norm sparsity penalty on the weight matrix. The mixed norm sparsity penalty is the glyph[lscript] 1 norm of the glyph[lscript] 2 norm on the weights corresponding to the same variable used over all the PLS components. The proposed global penalty encourages the selected variables to be shared among all the K PLS components. A novel augmented Lagrangian method is proposed to solve the optimization problem, which enables us to obtain the global SIMPLS components and to perform joint variable selection simultaneously. A greedy algorithm is proposed to overcome the computation difficulties in the iterations, and soft thresholding for sparsity occurs naturally as part of the iterative solution. Experiments show that our approach to PLS regression attains better or as good performance (lower mean squared error, MSE) with many fewer selected predictor variables. These experiments include a chemometric data set, and a human viral challenge study dataset, in addition to numerical simulations.

We review the developments in PLS-R for both univariate and multivariate responses in Section 2, in which we discuss different objective functions that have been proposed for PLS-R, particularly the Statistically Inspired Modification of PLS (SIMPLS) proposed by de Jong [10]. In Section 3, we formulate the Jointly Sparse Global SIMPLS-R by proposing a new criterion that jointly optimizes over K weight vectors (components) and imposing a mixed norm sparsity penalty to select variables jointly. The algorithmic implementation is discussed in Section 4 with simulation experiments presented in Section 5. The proposed Jointly Sparse Global SIMPLS Regression is applied to two applications: (1) Chemometrics in Section 6 and (2) Predictive health studies in Section 7. Section 8 concludes the paper.

## 2. Partial Least Squares Regression

Partial Least Squares (PLS) methods embrace a suite of data analysis techniques based on algorithms belonging to the PLS family. These algorithms consist of various extensions of the Nonlinear estimation by Iterative PArtial Least Squares (NIPALS) algorithm that was proposed by Herman Wold [11] as an alternative algorithm for implementing Principal Component Analysis (PCA) [12]. The NIPALS approach was slightly modified by Svante Wold, and Harald Martens, in order to obtain a regularized component based regression tool, known as PLS Regression (PLS-R) [1, 8].

Suppose that the data consists of n samples of p independent variables X ∈ R n × p and q dependent variables (responses) Y ∈ R n × q . In standard PLS Regression the aim is to define orthogonal latent components in R p , and then use

such latent components as predictors for Y in an ordinary least squares framework. The X weights used to compute the latent components can be specified by using iterative algorithms belonging to the NIPALS family or by a sequence of eigen-decompositions. The general underlying model is X = TP ′ + E and Y = TQ ′ + F , where T ∈ R n × K is the latent component matrix, P ∈ R p × K and Q ∈ R q × K are the loading matrices, K is the number of components, E and F are the residual terms. The latent components in T = [ t 1 t 2 ... t K ] are linear combinations of the independent variables, hence PLS can be viewed as a dimensional reduction technique, reducing the dimension from p to K . The latent components should be orthogonal to each other either by construction as in NIPALS [1, 8] or via constrained optimizations as in SIMPLS [10]. This allows PLS to build a parsimonious model for high dimensional data with collinearity [8].

## 2.1. Univariate Response

We assume, without loss of generality, that all the variables have been centered in a pre-processing step. For univariate Y , i.e q = 1, PLS Regression, also often denoted as PLS1, successively finds X weights R = [ r 1 r 2 ... r K ] as the solution to the constrained optimization

$$\mathbf r _ { k } = \arg \max _ { \mathbf r } \{ \mathbf r ^ { \prime } X ^ { \prime } _ { ( k - 1 ) } Y Y ^ { \prime } X _ { ( k - 1 ) } \mathbf r \} \ s. t. \ \mathbf r ^ { \prime } \mathbf r = 1. \quad \quad ( 2. 1 )$$

where X ( k -1) is the matrix of the residuals (i.e., the deflated matrix) from the regression of the X -variables on the first k -1 latent components, and X (0) = X . These weights are then used to find the latent components T = [ X (0) r 1 X (1) r 2 ... X ( K -1) r K ]. Such components can be also expressed in terms of original variables (instead of deflated variables), i.e., as T = XW , where W is the matrix containing the weights to be applied to the original variables in order to exactly obtain the latent components [13].

For a fixed number of components, the response variable Y is predicted in an ordinary least squares regression model, where the latent components play the role of the exogenous variables,

$$\hat { Q } = \arg \min _ { Q } \{ | | Y - T Q ^ { \prime } | | _ { 2 } \} = ( T ^ { \prime } T ) ^ { - 1 } T ^ { \prime } Y.$$

This provides the regression coefficients ˆ β PLS = WQ ˆ ′ for the model Y = Xβ PLS + F .

Depending on the number of selected latent components the length ‖ ˆ β PLS ‖ 2 of the vector of PLS coefficients changes. In particular, de Jong [14] had shown that the sequence of these coefficient vectors has lengths that are strictly increasing as the number of components increases. This sequence converges to the ordinary least squares coefficient vector and the maximum number of latent components obtainable equals the rank of the X matrix. Thus, by using a

number of latent components K &lt; rank X ( ), PLS-R performs a dimension reduction by shrinking the X matrix. Hence, PLS-R is a suitable tool for problems for which the data contains many more variables p than observations n .

The objective function in (2.1) can be interpreted as maximizing the squared covariance between Y and the latent component: corr 2 ( Y, X k -1 r k )var( X k -1 r k ). Because the response Y has been taken into account to formulate the latent matrix, PLS usually has better performance in prediction problems than principle component analysis (PCA) does. This is one of the main differences between PLS-R and PCA [15].

## 2.2. Multivariate Response

Similarly to univariate response PLS-R, multivariate response PLS-R selects latent components in R p and R q , i.e., t k and v k , such that the covariance between t k and v k is maximized. For a specific component, the sets of weights r k ∈ R p and c k ∈ R q are obtained by solving

$$& \max \{ \mathfrak { t } ^ { \prime } \mathbf v \} = \max \{ \mathfrak { r } ^ { \prime } X ^ { \prime } _ { ( k - 1 ) } Y _ { ( k - 1 ) } \mathfrak { c } \} \\ & s. t. \ \, \mathfrak { r } ^ { \prime } \mathbf r = \mathfrak { c } ^ { \prime } \mathbf c = 1$$

where t k = X ( k -1) r k , v k = Y ( k -1) c k , and X ( k -1) and Y ( k -1) are the deflated matrices associated with X and Y . Notice that the optimal solution c k should be proportional to Y ′ k -1 X k -1 r k . Therefore, the optimization in (2.2) is equivalent to

$$& \max _ { r } \left \{ { \mathbf r } ^ { \prime } X ^ { \prime } _ { k - 1 } Y _ { k - 1 } Y ^ { \prime } _ { k - 1 } X _ { k - 1 } { \mathbf r } \right \} \\ & s. t. \ r ^ { \prime } { \mathbf r } = 1.$$

For each component, the solution to this criterion can be obtained by using a so called PLS2 algorithm. A detailed description of the iterative algorithm as presented by H¨skuldsson [16] is in Algorithm 1 . o

In 1993 de Jong proposed a variant of the PLS2 algorithm, called Statistically Inspired Modification of PLS (SIMPLS), which calculates the PLS latent components directly as linear combinations of the original variables [10]. The SIMPLS was first developed as the solution to an optimization problem

$$& \text{w} _ { k } = \underset { w } { \arg \max } ( \text{w} ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X \text{w} ) & & ( 2. 4 ) \\ & s. t. \quad \text{w} ^ { \prime } \text{w} = 1, \quad \text{w} ^ { \prime } X ^ { \prime } X \text{w} _ { j } = 0 \ f o r \ j = 1, \dots, k - 1.$$

Ter Braak and de Jong [17] provided a detailed comparison between the objective functions for PLS2 in (2.3) and SIMPLS in (2.4) and showed that the successive weight vectors w k can be derived either from the deflated data matrices or the original variables in PLS2 and SIMPLS respectively. Let W + be the Moore-Penrose inverse of W = [ w w 1 2 ... w k -1 ]. The PLS2 algorithm (Algorithm 1) is equivalent to solving the optimization

$$w _ { k } = \underset { w } { \arg \max } ( w ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X w )$$

## Algorithm 1: PLS2 algorithm

```
  Algorithm  1: PLS2 algorithm
    
    for k=1:K do
      |    initialize r
        X = X.new
        Y = Y.new
        while solution has not converged do
          |    t = Xr
          |    c = Y't
          |    Scale c to length 1
          |    v = Yc
          |    r = X'v
          |    Scale r to length 1
          loading vector p = X't/(t't)
        deflate X.new = X - tp'
        regression b = Y't/(t't)
        deflate Y.new = Y - tb'
      |    r.k = r
```

$$s. t. \, { \mathbf w } ^ { \prime } ( I - W W ^ { + } ) { \mathbf w } = 1, { \mathbf w } ^ { \prime } X ^ { \prime } X { \mathbf w } _ { i } = 0 \ f o r \ i = 1, \dots, k - 1.$$

Both NIPALS and SIMPLS have the same objective function but each is maximized under a different normalization constraint. NIPALS and SIMPLS are equivalent when Y is univariate, but provide slightly different weight vectors in multivariate scenarios. The performance depends on the nature of the data, but SIMPLS appears easier to interpret since it does not involve deflation of the data set [10]. We develop our globally sparse PLS-R based on the SIMPLS optimization formulation.

## 3. Sparse Partial Least Squares Regression

## 3.1. glyph[lscript] 1 Penalized Sparse PLS Regression

One approach to sparse PLS-R is to add the glyph[lscript] 1 norm of the weight vector, a sparsity inducing penalty, to (2.4). The solution for the first component would be obtained by solving

$$& w _ { 1 } = \underset { w } { \arg \max } ( \mathbf w ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X \mathbf w ) & & ( 3. 1 ) \\ & s. t. \ \mathbf w ^ { \prime } \mathbf w = 1, \ | | \mathbf w | | _ { 1 } \leq \lambda.$$

The addition of the glyph[lscript] 1 norm is similar to SCOTLASS (simplified component lasso technique), the sparse PCA proposed by Jolliffe [18]. However, the solution of SCOTLASS is not sufficiently sparse, and the same issue remains in (3.1). Chun and Keles [9] reformulated the problem, promoting the exact zero property by imposing the glyph[lscript] 1 penalty on a surrogate of the weight vector instead of the original weight vector [9]. For the first component, they solve the following optimization by alternating between updating w and updating z (block

coordinate descent).

$$w _ { 1 }, z _ { 1 } = \underset { w, z } { \arg \min } \{ - \kappa w ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X w + ( 1 - \kappa ) ( z - w ) ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X ( z - w ) + \lambda _ { 1 } \| z \| _ { 1 } + \lambda _ { 2 } \| z \| _ { 2 } ^ { 2 } \}$$

$$s. t. \ w ^ { \prime } w = 1$$

Allen et. al proposed a general framework for regularized PLS-R [20].

$$\max _ { w, v } \, w ^ { \prime } M \mathbf v - \lambda P ( \mathbf w ) \ s. t. \ w ^ { \prime } \mathbf w \leq 1, \mathbf v ^ { \prime } \mathbf v = 1$$

in which M is the cross-product matrix X Y ′ , and the regularization function P is a convex penalty function. The formulation is a relaxation of SIMPLS with penalties being applied to the weight vectors, and can be viewed as a generalization of [9].

As mentioned in the Introduction, these formulations ([9, 20]) penalize the variables in each PLS component independently. This paper proposes an alternative in which variables are penalized simultaneously over all components. First, we define the global weight matrix, consisting of the K weight vectors, as

$$W = \left [ \begin{array} { c c c c } | & | & | & & | \\ w _ { 1 } & w _ { 2 } & \cdots & w _ { K } \\ | & | & | & & | \end{array} \right ] = \left [ \begin{array} { c c c$$

Notice that the elements in a particular row of W, i.e., w ′ ( j ) , are all associated with the same predictor variable x j . Therefore, rows of zeros correspond to variables that are not selected. To illustrate the drawbacks of penalizing each variable in each component independently, as in [9], suppose that each entry in W is selected independently with probability p 1 . The probability that the ( j ) th variable is not selected becomes (1 -p 1 ) K , and the probability that all the variables are selected by at least one weight vector is [1 -(1 -p 1 ) K p ] , which increases as the number of weight vectors K increases. This suggests that for large K the local variable selection approach of [9] may not lead to an overall sparse and parsimonious PLS-R model. In such cases a group sparsity constraint can be employed to limit the overall number of selected variables.

## 3.2. Jointly Sparse Global SIMPLS Regression

The Jointly Sparse Global SIMPLS Regression variable selection problem is to find the top K weight vectors that best relate X to Y , while using limited number of variables. This is a subset selection problem that is equivalent to adding a constraint on the glyph[lscript] 0 norm of the vector consisting of the norms of the rows of W , i.e, the number of nonzero rows in W . For concreteness we use the glyph[lscript] 2 norm for the rows. This leads to the optimization problem

$$W = \underset { W } { \arg \min } - \frac { 1 } { n ^ { 2 } } \sum _ { k = 1 } ^ { K } w _ { k } ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X w _ { k } \quad \quad \quad$$

glyph[negationslash]

s.t. || glyph[pi1] || 0 ≤ t, w w ′ k k = 1 ∀ k, and w ′ k X X ′ w i = 0 ∀ i = k in which

$$\varpi = \left [ \begin{array} { c } | | \mathbf w _ { ( 1 ) } | | _ { 2 } \\ | | \mathbf w _ { ( 2 ) } | | _ { 2 } \\ \vdots \\ | | \mathbf w _ { ( p ) } | | _ { 2 }$$

The objective function (3.2), which we refer to as global SIMPLS, is the sum of the objective functions (2.4) in the first K iterations of SIMPLS. Instead of the sequential greedy solution in PLS2 algorithm, the proposed jointly sparse global SIMPLS Regression solves for the K weight vectors simultaneously. We introduce the 1 n 2 factor to the objective function to interpret it as an empirical covariance. Given the complexity of this combinatorial problem, as in standard optimization practice, we relax the glyph[lscript] 0 norm optimization to a mixed norm structured sparsity penalty [21].

glyph[negationslash]

$$W = & \underset { W } { \arg } \min - \frac { 1 } {$$

The glyph[lscript] 2 norm of each row of W promotes grouping entries in W that relate to the same predictor variable, whereas the glyph[lscript] 1 norm promotes a small number of groups, as in (3.1).

## 4. Algorithmic Implementation for Jointly Sparse Global SIMPLS Regression

Constrained eigen-decomposition and group variable selection are each wellstudied problems for which efficient algorithms have been developed. We propose to solve the optimization (3.3) by augmented Lagrangian methods, which allows one to solve (3.3) by variable splitting iterations. Augmented Lagrangian methods introduce a new variable M , constrained such that M = W , such that the row vectors m ( j ) of M obey the same structural pattern as the rows of W :

glyph[negationslash]

$$& \min _ { W, M } - \frac { 1 } { n ^ { 2 } }$$

The optimization (4.1) can be solved by replacing the constrained problem by an unconstrained one with an additional penalty on the Frobenius norm of the

difference M -W . This penalized optimization can be iteratively solved by an alternating direction method of multipliers (ADMM) algorithm [25, 26, 27, 28, 29, 30, 31, 32], a block coordinate descent method that alternates between optimizing over W and over M (See Algorithm 2). We initialize Algorithm 2 with M (0) equal to the solution of SIMPLS, and D (0) equal to the zero matrix. Setting the parameter µ is nontrivial [31], and is hand-tuned for fastest convergence in some applications [27]. Once the algorithm converges, the final PLS regression coefficients are obtained by applying SIMPLS regression on the selected variables keeping the same number of components K . The optimization over W can be further simplified to a secular equation problem, whereas the optimization over M can be shown to reduce to a soft thresholding operation. The algorithm iterates until the stopping criterion based on the norm of the residuals || W ( τ ) -M ( τ ) || F &lt; glyph[epsilon1] is satisfied, for some given tolerance glyph[epsilon1] . As described later in the experimental comparisons section, the parameters λ and K are decided by cross validation.

Algorithm 2: Algorithm for solving the global SIMPLS with global variable selection problem using the augmented Lagrangian method.

glyph[negationslash]

$$\begin{array} {$$

Optimization over W The following optimization in Algorithm 2 is a nonconvex quadratically constrained quadratic program (QCQP).

$$W ^ { ( \tau + 1 ) } = \underset { W } { \arg \$$

glyph[negationslash]

$$s. t. \, { \mathbf w } ^ { \prime } _ { k } {$$

We propose solving for the K vectors in W successively by a greedy approach. Let m k and d k be the columns of the matrices M and D , and ω k = m k + d k . The optimization over W becomes

$$\mathbf w _ { k } ^ { ( \tau + 1 ) } \quad = \quad & \underset { \mathbf w } { \arg \min } - \frac { 1 } { n ^ { 2 } } \mathbf w ^ { \prime } X ^ { \prime } Y Y ^ { \prime } X \mathbf w + \frac { \mu } { 2 } \| \mathbf w - \omega _ { k } \| _ { 2 } ^ { 2 } } \\ s. t. \quad & \ \mathbf w ^ { \prime } \mathbf w = 1, \ \mathbf w ^ { \prime } X ^ { \prime } X \mathbf w _ { i } = 0 \ \forall \ i < k.$$

Lemma 4.1. Let N be an orthonormal basis for the orthogonal complement of { X X ′ w i } , i &lt; k . The optimization (4.2) can be solved by the method of Lagrange multipliers. The solution is w k = N A ( -αI ) -1 b , in which A = -1 n 2 N X Y Y XN ′ ′ ′ , b = µ 2 N ′ ω k and α is the minimum solution that satisfies b ′ ( A -αI ) -2 b = 1 .

Proof. Let w = N ˜ w , then the optimization (4.2) can be written as

$$\min \tilde { w } ^ { \prime } A \tilde { w } - 2 b ^ { \prime } \tilde { w } \ s. t. \ \tilde { w } ^ { \prime } \tilde { w } = 1.$$

Since we assume that w is a linear combination of the basis vectors in N , the orthogonality conditions in (4.2) are automatically satisfied. Hence these conditions have been dropped in the new formulation. Then using Lagrange multipliers, we can show that the solution takes the form as stated above.

Suppose there are two solutions of α that satisfy b ′ ( A -αI ) -2 b = 1, corresponding to two pairs of solutions to the optimization, ( ˜ w 1 , α 1 ) and ( ˜ w 2 , α 2 ). Since ˜ = ( w A -αI ) -1 b ,

$$A \tilde { w } _ { 1 } = \alpha _ { 1 } \tilde { w } _ { 1 } + b$$

$$A \tilde { w } _ { 2 } = \alpha _ { 2 } \tilde { w } _ { 2 } + b.$$

By multiplying (4.3) by ˜ w 1 , and (4.4) by ˜ w 2 , then subtracting the two new equations, we have

$$\tilde { w } _ { 1 } ^ { \prime } A \tilde { w } _ { 1 } - \tilde { w } _ { 2 } ^ { \prime } A \tilde { w } _ { 2 } = \alpha _ { 1 } - \alpha _ { 2 } + ( \tilde { w } _ { 1 } ^ { \prime } - \tilde { w } _ { 2 } ^ { \prime } ) b.$$

On the other hand, by multiplying (4.3) by ˜ w 2 , and (4.4) by ˜ w 1 , and subtracting the two new equations, we have

$$( \tilde { w } _ { 1 } ^ { \prime } - \tilde { w } _ { 2 } ^ { \prime } ) b = ( \alpha _ { 1 } - \alpha _ { 2 } ) \tilde { w } _ { 1 } ^ { \prime } \tilde { w } _ { 2 }.$$

Given (4.5) and (4.6), it can be shown that

$$( \tilde { w } _ { 1 } ^ { \prime } A \tilde { w } _ { 1 } - 2 b ^ { \prime } \tilde { w } _ { 1 } ) - ( \tilde { w } _ { 2 } ^ { \prime } A \tilde { w } _ { 2 } - 2 b ^ { \prime } \tilde { w } _ { 2 } ) = ( \alpha _ { 1 } - \alpha _ { 2 } ) - ( \tilde { w } _ { 1 } ^ { \prime } - \tilde { w } _ { 2 } ^ { \prime } ) b = \frac { \alpha _ { 1 } - \alpha _ { 2 } } { 2 } | | \tilde { w } _ { 1 } - \tilde { w } _ { 2 } | | _ { 2 } ^ { 2 }.$$

Hence, one should select the minimum among all the feasible α

's.

The equation b ′ ( A -αI ) -2 b = 1 is a secular equation, a well studied problem in constrained eigenvalue decomposition [23, 24]. The more general problem of least squares with a quadratic constraint was discussed in [22]. We can diagonalize the matrix A as A = UDU ′ , in which D is diagonal with eigenvalues d , d 1 2 , ..., d p in decreasing order on the diagonal, and the columns of U are the corresponding eigenvectors. Define

$$g ( \alpha ) = b ^ { \prime } ( A - \alpha I ) ^ { - 2 } \mathbf b = b ^ { \prime } U \left [ \begin{array} { c c c c } \frac { 1 } { ( d _ { 1 } - \alpha ) ^ { 2 } } & & & \\ & \frac { 1 } { ( d _ { 2 } - \alpha ) ^ { 2 } } & \\ & & \ddots & \\ & & & \frac { 1 } { ( d _ { p } - \alpha ) ^ { 2 } } \end{array} \right ] U ^ { \prime } \mathbf b.$$

Let ˜ b = U ′ b , then g α ( ) = b ′ ( A -αI ) -2 b = ∑ i ˜ b 2 i ( d i -α ) 2 , and hence g α ( ) = 1 is a secular equation. g α ( ) increases strictly as α increases from -∞ to d p , since

$$g ^ { \prime } ( \alpha ) = \sum _ { i } \frac { 2 \tilde { b } _ { i } ^ { 2 } } { ( d _ { i } - \alpha ) ^ { 3 } }$$

is positive for -∞ &lt; α &lt; d p . Moreover, given the limits

$$\lim _ { \mu \mu \mu \nu } g ( \alpha ) = 0$$

$$\lim _ { \alpha \to d _ { p } ^ { - } } g ( \alpha ) = \infty$$

$$\underset { \lambda - 1 } { \min } g ( \alpha ) = \upsilon \\ \lim _ { \lambda - 1 } g ( \alpha ) = \infty$$

we can conclude that there is exactly one solution α &lt; d p to the equation g α ( ) = 1, [23]. An iterative algorithm (Algorithm 3 ) is used to solve g α ( ) = 1 starting from a point to the left of the smallest eigenvalue d p [24].

Notice that calculating g α ( ) = b ′ ( A -αI ) -2 b involves inverting a ( p -k +1) × ( p -k +1) matrix, but A has rank at most q . We can reduce the computational burden by the use of Woodbury matrix identity,

$$( A - \alpha I ) ^ { - 1 } = \frac { - 1 } { \alpha } ( I - \frac { 1 } { \alpha n ^ { 2 } } N ^ { \prime } X ^ { \prime } Y ( I + \frac { 1 } { \alpha n ^ { 2 } } Y ^ { \prime } X N N ^ { \prime } X ^ { \prime } Y ) ^ { - 1 } Y ^ { \prime } X N ).$$

The new format only requires inverting a q × q matrix, and in most applications, the number of responses q is much less than the number of predictors p . Furthermore, N is involved in g α ( ) in the form of NN ′ = I -HH ′ , in which H is an orthonormal basis for { X X ′ w i } , i &lt; k . H can be constructed by Gram-Schmidt process as the algorithm successively finds the weight vectors w ′ i s .

## Algorithm 3: Iteration for solving secular equation.

$$\begin{matrix} \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \	 & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &\quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \null & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &
 & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \\ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &
 \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ } & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &  \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &
 \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \ \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \ \quad & \quad & \ \quad & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \quad & \quad & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \ & \quad & \quad & \quad & \ \quad & \ \quad & \ \quad & \quad & \ \quad & \ \quad & \quad & \quad & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \quad & \ \quad & \quad & \ \quad & \ \ } & \quad & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \ } & \quad & \ \quad & \ \ } & \ \ \ \\ & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \ } & \quad & \ \quad & \ \quad & \ \quad & \ \quad & \ \ \\ & \quad & \quad & \ \ \ & \quad & \ \ & \quad & \ \ \ & & \quad & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ & \ \ & & \ \ & \ & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & \ & \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ & & \ \ & \ \ & & \ \ & \ & \ \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & \ \ & \ & & \ \ & \ \ & \ \ & & \ \ & \ & \ \ & & \ \ & \ & \ \ & & \ \ & \ & \ & \ \ & \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & \ & \ \ & & \ \ & \ \ & \ \ & & \ \ & \ \ & \ & & \ \ & \ & \ & & \ \ & & \ \ & \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & & \ \ & & \ \ & \ \ & & \ \ & \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & \ & \ & \ & \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ & \ \ & & \ & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ \ & & \ & & \ & & \ & \ & & \ \$$

Optimization over M The optimization over M has a closed form solution. Let ∆ = W ( τ +1) -D ( τ ) , and δ ( j ) denote the j th row of ∆, then each row of M is given as m ( j ) = [ || δ ( j ) || -λ µ ] + δ ( j ) || δ ( j ) || , in which [ z ] + = max z, { 0 } .

Convergence analysis of ADMM can be found in [25, 29, 30, 32]. In particular, it has been shown that ADMM converges linearly for strongly convex objective functions [30]. Although the convergence is based on strong convexity assumptions, ADMM has been widely applied in practice [26, 27], even to nonconvex problems [28]. The proposed Jointly Sparse Global SIMPLS Regression is one example of these nonconvex applications.

## 5. Simulation Experiments

We implement the simulation models in [33, 9]. There are four models all following Y = Xβ + F , in which the number of observation is n = 100, and the dimension is p = 5000. Full details of the models are given in Table 1. We compare five different methods: the standard PLS regression (denoted as PLS in the following comparison tables), PLS generalized linear regression proposed by Bastien et al. [34], glyph[lscript] 1 penalized PLS regression [9] (denoted as glyph[lscript] 1 SPLS), Lasso [35] and the Jointly Sparse Global SIMPLS regression ( denoted as glyph[lscript] 1 /glyph[lscript] 2 SPLS). All the methods select their parameters by ten fold cross-validation, except for the PLS generalized linear regression, which stops including an additional component if the new component is not significant. The parameter µ in the Jointly Sparse Global SIMPLS-R is fixed to 2000, and updated in each iteration by a scaling factor 1 01 . µ . The experiments on real data in Section 6 and 7 are using the same setting of µ . Two i.i.d sets are generated for each trial: one as the training set and one as the test set. Ten trials are conducted for each model, and the averaged results are listed in Table 2.

In most of the simulations, we observe that the proposed Jointly Sparse Global SIMPLS-R performs as good or better than other methods in terms of the prediction MSE. In particular, the number of variables and the number of components chosen in Jointly Sparse Global SIMPLS-R are usually less than the glyph[lscript] 1 penalized PLS-R. We also calculate the R 2 for each method on the training data to measure the variation explained. The standard PLS regression and the PLS generalized linear regression proposed by Bastien et al. [34] both have R 2 close to 1, but the performance in terms of MSE is not ideal in the first three models for these methods. This may suggest that these models overfit the data. In addition to the averaged performance, the p-values of one sided paired t-test suggest that Jointly Sparse Global SIMPLS-R reduces model complexity significantly from those in the standard PLS regression and the PLS generalized linear regression. Lasso achieves low complexity in terms of the number of variables, but the MSE is high compared to Jointly Sparse Global SIMPLS-R and the glyph[lscript] 1 penalized PLS-R. The cross validation time for Jointly Sparse Global SIMPLS Regression is long, searching over a two-dimensional grid of the number of components K and the regularization parameter λ to minimize MSE. However, the performance in terms of prediction MSE improves, and the model complexity in terms of the number of variables and the number of components both decreases compared with other methods in most simulations.

## Table 1

Simulation models in [33, 9]: These models were originally proposed to test supervised principal components. Model 1 and Model 2 are designed such that one latent variable dominates the multi-collinearity. Model 3 has multiple latent variables, and Model 4 has a correlation structure from an autoregressive process. In these model descriptions, I is used to denote the indicator function. We use i ( 1 ≤ i ≤ n ) to index the n samples, j ( 1 ≤ j ≤ p ) to index the p variables, and k to index the hidden components.

## Model 1

H 1 ( ) = 3 i I ( i ≤ 50) + 4 I ( i &gt; 50) , 1 ≤ i ≤ n

H 2 ( ) = 3 5 i .

X j = H k + ε , j p k - 1 &lt; j ≤ p k , k = 1 2 , , ( p , ..., p 0 2 ) = (0 , 50 , p )

β

j

=

{

1

25

1

≤

j

≤

50

0

51

≤

j

≤

p

ε j is N (0 , I n ) distributed, and F is N (0 , 1 5 . 2 I n ) distributed.

## Model 2

H 1 ( ) = 2 5 i . I ( i ≤ 50) + 4 I ( i &gt; 50)

H 2 ( ) = 3 5 + 1 5 i . . I ( u 1 i ≤ 0 4) .

H 3 ( ) = 3 5 + 0 5 i . . I ( u 2 i ≤ 0 7) .

H 4 ( ) = 3 5 i . - 1 5 . I ( u 3 i ≤ 0 3) .

H 5 ( ) = 3 5 i .

u 1 i , u 2 i , u 3 i are i.i.d from Unif (0 , 1)

X j = H k + ε , j p k - 1 &lt; j ≤ p k , k = 1 , ..., 5 ( , p , ..., p 0 5 ) = (0 , 50 100 200 300 , , , , p )

β

j

=

{

1

25

1

≤

j

≤

50

0

51

≤

j

≤

p

ε j is N (0 , I n ) distributed, and F is N (0 , I n ) distributed.

## Model 3

H 1 ( ) = 2 5 i . I ( i ≤ 50) + 4 I ( i &gt; 50)

H 2 ( ) = 2 5 i . I (1 ≤ i ≤ 25 , or 51 ≤ i ≤ 75) + 4 I (26 ≤ i ≤ 50 , or 76 ≤ i ≤ 100)

H 3 ( ) = 3 5 + 1 5 i . . I ( u 1 i ≤ 0 4) .

H 4 ( ) = 3 5 + 0 5 i . . I ( u 2 i ≤ 0 7) .

H 5 ( ) = 3 5 i . - 1 5 . I ( u 3 i ≤ 0 3) .

H 6 ( ) = 3 5 i .

u 1 i , u 2 i , u 3 i are i.i.d from Unif (0 , 1)

X j = H k + ε , j p k - 1 &lt; j ≤ p k , k = 1 , ..., 5 ( , p , ..., p 0 6 ) = (0 , 25 50 100 200 300 , , , , , p )

β j = 25 0 51 ≤ j ≤ p

{

1

1

≤

j

≤

50

ε j is N (0 , I n ) distributed, and F is N (0 , I n ) distributed.

## Model 4

H 1 ( ) = i I ( i ≤ 50) + 6 I ( i &gt; 50)

H 2 ( ) = 3 5 + 1 5 i . . I ( u 1 i ≤ 0 4) .

H

3

( ) = 3 5 + 0 5

i

.

.

I

(

u

2

i

≤

0 7)

.

H 4 ( ) = 3 5 i . - 1 5 . I ( u 3 i ≤ 0 3) .

H 5 ( ) = 3 5 i .

u 1 i , u 2 i , u 3 i are i.i.d from Unif (0 , 1)

X = ( X (1) , X (2) )

X (1) is generated from N (0 , Σ 50 × 50 ) , Σ is from AR (1) with ρ = 0 9 . .

X (2) j = H k + ε , j p k - 1 &lt; i ≤ p k , k = 1 , ..., 5 ( , p , ..., p 0 5 ) = (0 , 50 100 200 300 , , , , p - 50)

β i = r m for p m - 1 &lt; i ≤ p m , m = 1 , ..., 6 , where

( p , ..., p 0 6 ) = (0 , 10 20 30 40 50 , , , , , p ) , ( r 1 , ..., r 6 ) = (8 , 6 4 2 1 0) , , , , / 25

ε i is N (0 , I n ) distributed, and F is N (0 , 1 5 . 2 I n ) distributed.

Table 2

Performance comparison table for the 4 simulation models. We compare five different methods: the standard PLS regression (PLS), PLS generalized linear regression proposed by Bastien et al. (Bastien), glyph[lscript] 1 penalized PLS regression ( glyph[lscript] 1 SPLS), Lasso and the Jointly Sparse Global SIMPLS regression ( glyph[lscript] 1 /glyph[lscript] 2 SPLS). These results are obtained by averaging over 10 trails.

| (5,4)                       | NA 0 . 053 0 . 006                                                                                   | NA 0 . 0119 2 . 70 × 10 - 4                                                                          | NA 0 . 1726 0 . 0054                                                                                 | NA 0 . 0485 0 . 0056                                                                                 |
|-----------------------------|------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------|
| paired t-test (5,3)         | 0 . 19 0 . 35 0 . 075                                                                                | 0 . 0184 0 . 1041 0 . 0344                                                                           | 0 . 3632 0 . 1697 0 . 3430                                                                           | 0 . 2201 0 . 3918 0 . 3812                                                                           |
| p-values of one sided (5,2) | 4 . 49 × 10 - 7 5 . 46 × 10 - 5 0 . 043                                                              | 1 . 19 × 10 - 11 6 . 91 × 10 - 9 0 . 0046                                                            | 5 . 53 × 10 - 6 1 . 39 × 10 - 12 0 . 0104                                                            | 6 . 82 × 10 - 5 0 . 4618 0 . 1874                                                                    |
| (5,1)                       | 0 . 5 1 . 84 × 10 - 11 0 . 01                                                                        | 0 . 0671 1 . 13 × 10 - 14 8 . 69 × 10 - 4                                                            | 0 . 4201 2 . 32 × 10 - 19 1 . 18 × 10 - 4                                                            | 0 . 4057 9 . 62 × 10 - 5 0 . 0087                                                                    |
| 1 /glyph[lscript] 2 SPLS    | 1.4 276.1 2.82 0.83 11414 6.40 0.002 11420                                                           | 1.1 171.7 2.69 0.75 11420 5.72 0.001 11426                                                           | 1.5 60.5 1.25 0.73 11295 5.48 0.001 11300                                                            | 2.1 1180.5 2.36 0.91 12051 7.97 0.004 12059                                                          |
| 4. Lasso                    | NA 40.7 3.23 0.59 51.59 0.04 0.04 51.67                                                              | NA 15.8 3.09 0.39 53.97 0.04 0.039 54.05                                                             | NA 41.3 1.48 0.75 49.45 0.04 0.039 49.53                                                             | NA 9.4 2.14 0.19 50.52 0.04 0.042 50.60                                                              |
| 3. glyph[lscript] 1 SPLS    | 1.9 246.5 3.00 0.71 43.49 0.05 0.002 43.54                                                           | 2.3 273.4 2.93 0.79 43.03 0.06 0.002 43.09                                                           | 1.4 89.2 1.27 0.77 43.84 0.04 0.001 43.88                                                            | 2.6 1260.8 2.41 0.78 44.31 0.10 0.007 44.41                                                          |
| 2. Bastien                  | 5 1129.4 2.98 1 0 121.91 0.011 121.92                                                                | 5 1158.4 2.99 1 0 122.69 0.011 122.70                                                                | 5 1156.4 1.48 1 0 126.08 0.01 126.09                                                                 | 5 1118.8 2.29 1 0 123.79 0.011 123.8                                                                 |
| 1. PLS                      | 1.4 5000 3.14 0.98 101.47 0.89 0.010 102.37                                                          | 2 5000 3.18 0.98 100.51 1.28 0.010 101.80                                                            | 1.4 5000 1.82 0.98 102.61 1.03 0.01 103.65                                                           | 2 5000 2.15 1 98.16 1.55 0.010 99.73                                                                 |
|                             | Model 1 number of comp. number of variables MSE R 2 Time CV Time analysis Time prediction Total time | Model 2 number of comp. number of variables MSE R 2 Time CV Time analysis Time prediction Total time | Model 3 number of comp. number of variables MSE R 2 Time CV Time analysis Time prediction Total time | Model 4 number of comp. number of variables MSE R 2 Time CV Time analysis Time prediction Total time |

Table 3 Performance comparison table for the Octane data.

|                                |         |                          |                                            | p-values of one sided paired t-test           | p-values of one sided paired t-test   |
|--------------------------------|---------|--------------------------|--------------------------------------------|-----------------------------------------------|---------------------------------------|
|                                | 1. PLS  | 2. glyph[lscript] 1 SPLS | 3. glyph[lscript] 1 /glyph[lscript] 2 SPLS | (3,1)                                         | (3,2)                                 |
| number of comp. number of var. | 5.5 225 | 4.5 87.3                 | 3.8 38.5                                   | 5 . 1888 × 10 - 16 1 . 1176 × 10 - 121 0.0032 | 0.0027 1 . 0967 × 10 - 15 0.1575      |
| MSE                            | 0.0564  | 0.0509                   | 0.0481                                     |                                               |                                       |

## 6. Application 1: Chemometrics Study

In this section we show experimental results obtained by comparing standard PLS-R, glyph[lscript] 1 penalized PLS-R [9] (denoted as glyph[lscript] 1 SPLS in the performance table), and our proposed Jointly Sparse Global SIMPLS-R (denoted as glyph[lscript] 1 /glyph[lscript] 2 SPLS in the performance table). All the methods have been applied on the Octane data set (see [13]). The Octane data is a real data set consisting of 39 gasoline samples for which the digitized Octane spectra have been recorded at 225 wavelengths (in nm). The aim is to predict the Octane number, a key measurement of the physical properties of gasoline, using the spectra as predictors. This is of major interest in real applications, because the conventional procedure to calculate the Octane number is time consuming and involves expensive and maintenanceintensive equipment as well as skilled labor.

The experiments are composed of 150 trials. In each trial we randomly split the 39 samples into 26 training samples and 13 test samples. The regularization parameter λ and number of components K are selected by 2-fold cross validation on the training set. The averaged results over the 150 trials are shown in Table 3. We further show the variable selection frequencies for the sparse PLS methods over the 150 trials superimposed on the octane data in Fig. 1 (B) and (C). In chemometrics, the rule of thumb is to look for variables that have large amplitudes in first derivatives with respect to wavelength. Notice that both glyph[lscript] 1 penalized PLS-R and Jointly Sparse Global SIMPLS-R have selected variables around 1200 and 1350 nm, and the selected region in the latter case is more confined. Box and Whisker plots for comparing the MSE, number of selected variables, and number of components of these three PLS formulations are shown in Fig. 1 (A). Comparing our proposed Jointly Sparse Global SIMPLS Regression with standard PLS-R and glyph[lscript] 1 penalized PLS-R [9], we show that Jointly Sparse Global SIMPLS-R attains better performance in terms of MSE, the number of predictors, and the number of components. Besides, the model complexity in Jointly Sparse Global SIMPLS-R is significantly lower than both standard PLS-R and glyph[lscript] 1 penalized PLS-R, given the p-values of one sided paired t-test.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0318v1-figure-2f891694-1264-4201-a52e-2a07504489b3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0318v1-figure-6b527d53-4675-4b8f-a184-d63bfb58217e.jpg)

(C)

Fig 1: (A) Box and Whisker plots for comparing the MSE, number of selected variables, and number of components of three PLS formulations: standard PLSR, glyph[lscript] 1 penalized PLS-R ( glyph[lscript] 1 SPLS), our proposed Jointly Sparse Global SIMPLSR ( glyph[lscript] 1 /glyph[lscript] 2 SPLS). (B) Variable selection frequency of glyph[lscript] 1 SPLS superimposed on the octane data and its first derivative: The height of the surfaces represents the exact value of the data over 225 variables for the 39 samples. The color of the surface shows the selection frequency of the variables as depicted on the colorbar. (C)Variable selection frequency of glyph[lscript] 1 /glyph[lscript] 2 SPLS superimposed on the octane data and its first derivative.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0318v1-figure-daea3dbe-7307-4e6f-9d51-86288f6510ed.jpg)

## 7. Application 2: Sparse Prediction of Disease Symptoms from Gene Expression

In this section we apply the Jointly Sparse Global SIMPLS Regression to 4 types of predictive health challenge studies involving the H3N2, the H1N1, the HRV, and the RSV viruses. In these challenge studies, publicly available from the NCBI-GEO website, serial peripheral blood samples were acquired from a population of subjects inoculated with live flu viruses [38, 39, 40, 41]. The prediction task in these experiments is to predict the symptom scores based on gene expression of 12023 genes. There were 10 symptom scores, i.e., runny nose, stuffy nose, sneezing, sore throat, earache, malaise, cough, shortness of breath, headache, and myalgia, documented over time. The symptoms are self-reported scores, ranging from 0 to 3. We linearly interpolate the gene expressions to match them with the sampling time of the symptom reports. We compare the Jointly Sparse Global SIMPLS-R with standard PLS-R and glyph[lscript] 1 penalized PLS-R by leaving one subject out as the test set, and the rest as the training set. The process is repeated until all subjects have been treated as the test set. The number of components for all methods and the regularization parameter in Jointly Sparse Global SIMPLS-R are selected by 2-fold cross validation to minimize the sum of the MSE of the responses. Since each subject has multiple samples, we perform the cross validation by splitting by subjects, i.e., no samples from the same subject will appear in both training and tuning sets. We restrict the responses to the first 3 symptoms, which are the upper respiratory symptoms, and the results are shown in Table 4. In most of the cases, the proposed Jointly Sparse Global SIMPLS-R method outperforms the standard PLS-R and glyph[lscript] 1 penalized PLS-R in terms of prediction MSE, number of components, number of genes. As can be seen in Table 4, the number of selected variables decreases significantly by applying the glyph[lscript] 1 /glyph[lscript] 2 mixed norm sparsity penalty to the PLSR objective function. Thus the proposed PLS-R method is able to construct a more parsimonious predictor relative to the other PLS-R methods having similar accuracy.

The PLS-R method can also be viewed as an exploratory data analysis tool for constructing low dimensional descriptors of the independent variables and response variables. Specifically, the general underlying matrix factorization model X = TP ′ + E and Y = TQ ′ + F , with latent component T = XW , provides a factor analysis model for the independent and response variables X and Y . T , P and Q can be interpreted in a similar manner as the singular vectors of PCA. However, different from PCA that does not account for the response variables, T , P and Q contain information about both the independent variables and the response variables. The factor analysis interpretation of the underlying PLS model is that T is a latent score matrix and P , Q are latent factor loading matrices that associate T with the independent variables and the response variables via the approximate matrix factorizations X ≈ TP ′ and Y ≈ TQ ′ , respectively. The correlations between the latent component T and the sum of the 3 upper respiratory symptoms are reported in Table 5, which also shows results for classic matrix factorization methods including non-negative matrix

Table 4

Performance comparison table for the predictive health study. We apply the Jointly Sparse Global SIMPLS Regression ( glyph[lscript] 1 /glyph[lscript] 2 SPLS) for sparse prediction of disease symptoms from gene expression. The performance is compared with standard PLS and glyph[lscript] 1 sparse PLS. glyph[lscript] 1 /glyph[lscript] 2 SPLS achieves lower MSE with significantly fewer variables in most of the studies.

|                      |        |                          |                                            | p-values of one sided paired t-test   | p-values of one sided paired t-test   |
|----------------------|--------|--------------------------|--------------------------------------------|---------------------------------------|---------------------------------------|
|                      | 1. PLS | 2. glyph[lscript] 1 SPLS | 3. glyph[lscript] 1 /glyph[lscript] 2 SPLS | (3,1)                                 | (3,2)                                 |
| H1N1 number of comp. | 2.8    | 2.3                      | 2.4                                        | 0.1163                                | 0.3322                                |
| number of genes      | 12023  | 3624.1                   | 3575.8                                     | 1 . 4451 × 10 - 10                    | 0.4842                                |
| Overall MSE          | 0.599  | 0.603                    | 0.591                                      | 0.2890                                | 0.1094                                |
| Runny nose MSE       | 0.167  | 0.165                    | 0.167                                      |                                       |                                       |
| Stuffy nose MSE      | 0.281  | 0.282                    | 0.269                                      |                                       |                                       |
| Sneezing MSE         | 0.151  | 0.157                    | 0.155                                      |                                       |                                       |
| H3N2 number of comp. | 3.2    | 2.5                      | 1.9                                        | 0.0030                                | 0.0863                                |
| number of genes      | 12023  | 3944.5                   | 1721.5                                     | 1 . 7547 × 10 - 9                     | 0.0601                                |
| Overall MSE          | 0.623  | 0.622                    | 0.609                                      | 0.3073                                | 0.2530                                |
| Runny nose MSE       | 0.186  | 0.174                    | 0.173                                      |                                       |                                       |
| Stuffy nose MSE      | 0.277  | 0.284                    | 0.272                                      |                                       |                                       |
| Sneezing MSE         | 0.160  | 0.164                    | 0.165                                      |                                       |                                       |
| HRV number of comp.  | 2.8    | 2.3                      | 2.2                                        | 0.0484                                | 0.3773                                |
| number of genes      | 12023  | 2193.2                   | 1779.1                                     | 5 . 5038 × 10 - 13                    | 0.3522                                |
| Overall MSE          | 0.628  | 0.607                    | 0.603                                      | 0.2020                                | 0.4490                                |
| Runny nose MSE       | 0.243  | 0.226                    | 0.232                                      |                                       |                                       |
| Stuffy nose MSE      | 0.324  | 0.323                    | 0.314                                      |                                       |                                       |
| Sneezing MSE         | 0.062  | 0.058                    | 0.057                                      |                                       |                                       |
| RSV number of comp.  | 3.2    | 2.3                      | 2.4                                        | 0.0198                                | 0.4103                                |
| number of genes      | 12023  | 2445.4                   | 3889.8                                     | 1 . 1584 × 10 - 9                     | 0.1472                                |
| Overall MSE          | 0.866  | 0.920                    | 0.855                                      | 0.3318                                | 0.0567                                |
| Runny nose MSE       | 0.312  | 0.327                    | 0.312                                      |                                       |                                       |
| Stuffy nose MSE      | 0.412  | 0.448                    | 0.397                                      |                                       |                                       |
| Sneezing MSE         | 0.143  | 0.145                    | 0.145                                      |                                       |                                       |

factorization (NMF) [42] and Bayesian linear unmixing (BLU) [43, 44], previously applied to this dataset, for comparison. Notice the sparse PLS-R methods achieve higher correlation, as expected. Remarkably, the proposed Jointly Sparse global SIMPLS-R achieves this higher degree of correlation with many fewer components and variables than the NMF and BLU methods. This experiment demonstrates that Jointly Sparse Global SIMPLS-R can be used as a factor analysis method to find the hidden molecular factors that best relate to the response.

## 8. Conclusion

The formulation of the global SIMPLS objective function with an added group sparsity penalty greatly reduces the number of variables used to predict the response. This suggests that when multiple components are desired, the variable selection technique should take into account the sparsity structure for the same variables among all the components. Our proposed Jointly Sparse Global SIMPLS Regression algorithm is able to achieve as good or better performance with fewer predictor variables and fewer components as compared to competing methods. It is thus useful for performing dimension reduction and variable

Table 5

Matrix factorization. We use the cross-validated parameters reported in Table 4 for standard PLS-R and glyph[lscript] 1 penalized PLS-R, and the Jointly Sparse Global SIMPLS-R to decide the number of components, and search over the grid { 1 2 , , ..., 10 } to find the number of factors that achieves the highest correlation for NMF and BLU. The correlation between each factor and the sum of responses is listed for each method. The first 3 methods, PLS, glyph[lscript] 1 SPLS, and glyph[lscript] 1 /glyph[lscript] 2 SPLS are supervised matrix factorizations, where as NMF and BLU are unsupervised. The unsupervised methods require many more factors to achieve comparable correlation.

| factor                                  | correlation of 1 2   |      | each factor 3   | with 4   | the sum 5   | of upper 6   | respiratory 7 8   |      | symptoms 9   |
|-----------------------------------------|----------------------|------|-----------------|----------|-------------|--------------|-------------------|------|--------------|
| H1N1                                    |                      |      |                 |          |             |              |                   |      |              |
| PLS                                     | 0.47                 | 0.38 | 0.38            |          |             |              |                   |      |              |
| glyph[lscript] 1 SPLS                   | 0.52                 | 0.37 |                 |          |             |              |                   |      |              |
| glyph[lscript] 1 /glyph[lscript] 2 SPLS | 0.57                 | 0.35 |                 |          |             |              |                   |      |              |
| NMF                                     | 0.32                 | 0.45 | 0.04            |          |             |              |                   |      |              |
| BLU                                     | 0.27                 | 0.19 | 0.02            | 0.51     | 0.16        | 0.06         |                   |      |              |
| H3N2                                    |                      |      |                 |          |             |              |                   |      |              |
| PLS                                     | 0.67                 | 0.42 | 0.33            |          |             |              |                   |      |              |
| glyph[lscript] 1 SPLS                   | 0.73                 | 0.33 | 0.33            |          |             |              |                   |      |              |
| glyph[lscript] 1 /L 2 SPLS              | 0.71                 | 0.33 |                 |          |             |              |                   |      |              |
| NMF                                     | 0.62                 | 0.70 | 0.10            |          |             |              |                   |      |              |
| BLU                                     | 0.54                 | 0.73 | 0.26            | 0.33     | 0.01        | 0.02         | 0.28              | 0.05 | 0.00         |
| HRV                                     |                      |      |                 |          |             |              |                   |      |              |
| PLS                                     | 0.45                 | 0.43 | 0.35            |          |             |              |                   |      |              |
| glyph[lscript] 1 SPLS                   | 0.52                 | 0.38 |                 |          |             |              |                   |      |              |
| glyph[lscript] 1 /glyph[lscript] 2 SPLS | 0.53                 | 0.42 |                 |          |             |              |                   |      |              |
| NMF                                     | 0.02                 | 0.22 | 0.19            |          |             |              |                   |      |              |
| BLU                                     | 0.11                 | 0.04 | 0.18            | 0.33     | 0.01        | 0.02         | 0.28              | 0.05 | 0.00         |
| RSV                                     |                      |      |                 |          |             |              |                   |      |              |
| PLS                                     | 0.66                 | 0.35 | 0.35            |          |             |              |                   |      |              |
| glyph[lscript] 1 SPLS                   | 0.70                 | 0.34 |                 |          |             |              |                   |      |              |
| glyph[lscript] 1 /glyph[lscript] 2 SPLS | 0.69                 | 0.39 |                 |          |             |              |                   |      |              |
| NMF                                     | 0.41                 | 0.13 | 0.16            | 0.01     | 0.31        | 0.02         | 0.11              | 0.11 | 0.01         |
| BLU                                     | 0.01                 | 0.02 | 0.23            | 0.20     | 0.03        | 0.68         | 0.12              |      |              |

selection simultaneously in applications with large dimensional data but comparatively few samples ( n &lt; p ).

The Jointly Sparse Global SIMPLS Regression objective function is minimized using augmented Lagrangian techniques and, in particular, the ADMM algorithm. The ADMM algorithm splits the optimization into an eigen-decomposition problem and a soft-thresholding that enforces sparsity constraints. The general framework is extendable to more complicated regularization and can thus be tailored for other PLS-type applications, e.g., positivity constraints or smoothness penalties. For example, in the chemometric application, the data is smooth over the wavelengths and we can apply wavelet shrinkage on the data or include a total variation regularization to encourage smoothness. The sparsity constraints can be imposed on the wavelet coefficients if wavelet shrinkage is applied, or together with total variation regularization. The equivalence of soft wavelet shrinkage and total variation regularization was discussed in [45]. One can also consider imposing sparsity structures on the weights corresponding to the same components, adding glyph[lscript] 1 penalty within the groups, or total variation regularization, depending on the applications. The decoupling property of the ADMM algorithm allows one to extend the Jointly Sparse Global SIMPLS Regression to these various regularizations.

## References

- [1] Wold, S., Martens, H., and Wold, H. (1983). The multivariate calibration problem in chemistry solved by the PLS method. Proceedings of the Conference on Matrix Pencils. Lectures Notes in Mathematics , 286-293.
- [2] Sj¨str¨m, M., and Wold, S., and Lindberg, W., and Persson, J., and Martens, o o H. (1983). A multivariate calibration problem in analytical chemistry solved by partial least-squares models in latent variables. Analytica Chimica Acta 150 , 61-70.
- [3] Martens, H., and Martens, M. (1999). Validation of PLS Regression models in sensory science by extended cross-validation. PLS'99 .
- [4] Rossouw, D., Robert-Grani´, C., and Besse, P. (2008). A sparse PLS for varie able selection when integrating omics data. Genetics and Molecular Biology , 7 (1), 35.
- [5] Chun, H., and Kele¸ s, S. (2009). Expression quantitative trait loci mapping with multivariate sparse partial least squares regression. Genetics , 182 (1), 79-90.
- [6] Chung, D., and Kele¸ s, S. (2010). Sparse partial least squares classification for high dimensional data. Statistical applications in genetics and molecular biology , 9 (1), 17.
- [7] Chun, H., Ballard, D. H., Cho, J., and Zhao, H. (2011). Identification of association between disease and multiple markers via sparse partial least squares regression. Genetic epidemiology , 35 (6), 479-486.
- [8] Wold, S., Ruhe, A., Wold, H., and Dunn III, W. J. (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.
- [9] Chun, H., and Kele¸ s, S. (2010). Sparse partial least squares regression for simultaneous dimension reduction and variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 72 , 3-25.
- [10] de Jong, S. (1993). SIMPLS: an alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems , 18 (3), 251263.
- [11] Wold, H. (1966). Nonlinear estimation by iterative least squares procedures. Research papers in statistics , 411-444.
- [12] Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of educational psychology , 24 (6), 417.
- [13] Tenenhaus, M. (1998). La R´ egression PLS: th´ eorie et pratique . Editions Technip.
- [14] de Jong, S. (1995). PLS shrinks. Journal of Chemometrics , 9 (4), 323-326.
- [15] Boulesteix, A. L., and Strimmer, K. (2007). Partial least squares: a versatile tool for the analysis of high-dimensional genomic data. Briefings in bioinformatics , 8 (1), 32-44.
- [16] H¨ oskuldsson, A. (1988). PLS regression methods. Journal of Chemometrics , 2 (3), 211-228.
- [17] ter Braak, C. J., and de Jong, S. (1998). The objective function of partial least squares regression. Journal of chemometrics , 12 (1), 41-54.
- [18] Jolliffe, I. T., Trendafilov, N. T., and Uddin, M. (2003). A modified principal component technique based on the LASSO. Journal of Computational and Graphical Statistics , 12 (3), 531-547.
- [19] Wold, H. (1975). Soft modelling by latent variables: the non-linear iterative partial least squares (NIPALS) approach. Perspectives in Probability and Statistics, In Honor of MS Bartlett , 117-144.
- [20] Allen, G. I., Peterson, C., Vannucci, M., and Maleti´-Savati´, M. (2013). c c Regularized partial least squares with an application to NMR spectroscopy. Statistical Analysis and Data Mining , 6 (4), 302-314.
- [21] Bach, F. R. (2008). Consistency of the group Lasso and multiple kernel learning. The Journal of Machine Learning Research , 9 , 1179-1225.
- [22] Gander, W. (1980). Least squares with a quadratic constraint. Numerische Mathematik , 36 (3), 291-307.
- [23] Gander, W., Golub, G. H., and von Matt, U. (1989). A constrained eigenvalue problem. Linear Algebra and its applications , 114 , 815-839.
- [24] Beck, A., Ben-Tal, A., and Teboulle, M. (2006). Finding a global optimal solution for a quadratically constrained fractional quadratic problem with applications to the regularized total least squares. SIAM Journal on Matrix Analysis and Applications , 28 (2), 425-445.
- [25] Eckstein, J., and Bertsekas, D. P. (1992). On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming , 55 (1-3), 293-318.
- [26] Goldstein, T., and Osher, S. (2009). The split Bregman method for L1regularized problems. SIAM Journal on Imaging Sciences , 2 (2), 323-343.

- [27] Afonso, M. V., Bioucas-Dias, J. M., and Figueiredo, M. A. (2011). An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems. Image Processing, IEEE Transactions on , 20 (3), 681-695.
- [28] Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends R © in Machine Learning , 3 (1), 1-122.
- [29] Hong, M., and Luo, Z. Q. (2012). On the linear convergence of the alternating direction method of multipliers. arXiv preprint arXiv:1208.3922.
- [30] Goldstein, T., O'Donoghue, B., and Setzer, S. (2012). Fast alternating direction optimization methods. CAM report, 12-35.
- [31] Ramani, S., and Fessler, J. A. (2012) A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction. Medical Imaging, IEEE Transactions on , 31 (3), 677-688.
- [32] Nien, H., and Fessler, J. A. (2014). A convergence proof of the split Bregman method for regularized least-squares problems. arXiv preprint arXiv:1402.4371.
- [33] Bair, E., Hastie, T., Paul, D., and Tibshirani, R. (2006). Prediction by supervised principal components. Journal of the American Statistical Association , 101 (473).
- [34] Bastien, P., Vinzi, V. E., and Tenenhaus, M. (2005). PLS generalised linear regression. Computational Statistics &amp; Data Analysis , 48 (1), 17-46.
- [35] Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological) , 267-288.
- [36] Martens, H., and Naes, T. (1989). Multivariate calibration . Wiley.
- [37] Masoum, S., Bouveresse, D. J. R., Vercauteren, J., Jalali-Heravi, M., and Rutledge, D. N. (2006). Discrimination of wines based on 2D NMR spectra using learning vector quantization neural networks and partial least squares discriminant analysis. Analytica chimica acta , 558 (1), 144-149.
- [38] Zaas, A. K., Chen, M., Varkey, J., Veldman, T., Hero, A. O., Lucas, J., Huang, Y., Turner, R., Gilbert, A., Lambkin-Williams, R., Øien, N. C., Nicholson, B., Kingsmore, S., Carin, L., Woods, C. W., and Ginsburg, G. S. (2009). Gene Expression Signatures Diagnose Influenza and Other Symptomatic Respiratory Viral Infections in Humans. Cell Host and Microbe , 6 (3), 207-217.
- [39] Huang, Y., Zaas, A. K., Rao, A., Dobigeon, N., Woolf, P. J., Veldman, T., Øien, N. C., McClain, M. T., Varkey, J. B., Nicholson, B., Carin, L., Kingsmore, S., Woods, C. W., Ginsburg, G. S., Hero, A. O. (2011). Temporal dynamics of host molecular responses differentiate symptomatic and asymptomatic influenza a infection. PLoS genetics , 7 (8), e1002234.
- [40] Woods., C. W., McClain, M. T., Chen, M., Zaas, A. K., Nicholson, B. P., Varkey, J., Veldman, T., Kingsmore, S. F., Huang, Y., Lambkin-Williams, R., Gilbert, A. G., Hero, A. O., Ramsburg, E., Glickman, S., Lucas1, J. E., Carin, L., and Ginsburg, G. S. (2013). A host transcriptional signature for presymptomatic detection of infection in humans exposed to influenza H1N1

or H3N2. PloS one 8.1: e52198.

- [41] Zaas, A. K., Burke, T., Chen, M., McClain, M., Nicholson, B., Veldman, T., Tsalik, E. L., Fowler, V., Rivers, E. P., Otero, R., Kingsmore, S. F., Voora, D., Lucas, J., Hero, A. O., Carin, L., Woods, C. W., and Ginsburg, G. S. (2013). A Host-Based RT-PCR Gene Expression Signature to Identify Acute Respiratory Viral Infection. Science translational medicine , 5 , 203ra126.
- [42] Paatero, P., Tapper, U. (1994). Positive matrix factorization: A nonnegative factor model with optimal utilization of error estimates of data values. Environmetrics , 5 (2), 111-126.
- [43] Dobigeon, N., Moussaoui, S., Coulon, M., Tourneret, J. Y., and Hero, A. O. (2009). Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery. Signal Processing, IEEE Transactions on , 57 (11), 4355-4368.
- [44] Bazot, C., Dobigeon, N., Tourneret, J.-Y., Zaas, A. K., Ginsburg, G. S. , Hero, A. O. (2013). Unsupervised Bayesian linear unmixing of gene expression microarrays. BMC Bioinformatics , 14 (1), 99.
- [45] Steidl, G., Weickert, J., Brox, T., Mrzek, P., and Welk, M. (2004). On the equivalence of soft wavelet shrinkage, total variation diffusion, total variation regularization, and SIDEs. SIAM Journal on Numerical Analysis , 42 (2), 686-713.
Tzu-Yu Liu Laura Trinchera Arthur Tenenhaus Dennis Wei Alfred O. Hero
Jointly Sparse Global SIMPLS Regression
Partial least squares (PLS) regression combines dimensionality reduction and
prediction using a latent variable model. Since partial least squares
regression (PLS-R) does not require matrix inversion or diagonalization, it can
be applied to problems with large numbers of variables. As predictor dimension
increases, variable selection becomes essential to avoid over-fitting, to
provide more accurate predictors and to yield more interpretable parameters. We
propose a global variable selection approach that penalizes the total number of
variables across all PLS components. Put another way, the proposed global
penalty encourages the selected variables to be shared among the PLS
components. We formulate PLS-R with joint sparsity as a variational
optimization problem with objective function equal to a novel global SIMPLS
criterion plus a mixed norm sparsity penalty on the weight matrix. The mixed
norm sparsity penalty is the $\ell_1$ norm of the $\ell_2$ norm on the weights
corresponding to the same variable used over all the PLS components. A novel
augmented Lagrangian method is proposed to solve the optimization problem and
soft thresholding for sparsity occurs naturally as part of the iterative
solution. Experiments show that the modified PLS-R attains better or as good
performance with many fewer selected predictor variables.