Weighted Sparse Partial Least Squares for Joint Sample and Feature Selection

Assume that there are two data matrices $\bm{X}\in\mathbb{R}^{n\times p}$ and $\bm{Y}\in\mathbb{R}^{n\times q}$ with a same set of samples and their columns are standardized to have mean zero and variance unit. PLS is equivalent to solving the following optimization problem [29]:

The basic idea of PLS is to obtain two projected vectors $\bm{u}$ and $\bm{v}$ by maximizing the covariance between the latent vectors $\bm{Xu}$ and $\bm{Yv}$. Specifically, a sparse PLS (sPLS) with $\ell_{0}$-norm constraint is formulated as:

where $k_{u}$ and $k_{v}$ are two parameters to control the sparsity of $\bm{u}$ and $\bm{v}$. The above model is also called sparse diagonal CCA problem [30, 31].

We observe that the objective function in Eq. (2) can be written as $\sum_{i=1}^{n}(\bm{Xu})_{i}(\bm{Yv})_{i}$. To consider the difference of samples, we set a weighted coefficient for each sample and modify the above objective function as $\sum_{i=1}^{n}w_{i}(\bm{Xu})_{i}(\bm{Yv})_{i}=\bm{u}^{T}\bm{X}^{T}\mbox{diag}(\bm{w})\bm{Y}\bm{v}$ and propose a weighted sparse PLS model with the $\ell_{\infty}/\ell_{0}$-norm constraints, which is denoted as ($\ell_{\infty}/\ell_{0}$-wsPLS or wsPLS for short):

where $\|\bm{w}\|_{\infty}=\max_{j}|w_{j}|$ and $\mbox{diag}(\bm{w})\in\mathbb{R}_{+}^{n\times n}$ is a diagonal matrix with the diagonal entries $w_{1},w_{2},\cdots,w_{n}$ indicating the weighted coefficients of samples and the $\ell_{\infty}/\ell_{0}$-norm constraints are used to select a subset of samples. The mathematical notations are summarized in Table I.

Eq. (3) is equivalent to the following optimization problem because $\bm{u}^{T}\bm{X}^{T}\mbox{diag}(\bm{w})\bm{Y}\bm{v}=\bm{w}^{T}[(\bm{X}\bm{u})\odot(\bm{Y}\bm{v})]$:

Obviously, the above $\ell_{\infty}/\ell_{0}$-wsPLS model reduces to sPLS model when all elements of $\bm{w}$ are one. We observe that $f(\bm{u},\bm{v},\bm{w})$ is also equal to $-\bm{u}^{T}\bm{X}^{T}[\bm{w}\odot(\bm{Y}\bm{v})]$ or $-\bm{v}^{T}\bm{Y}^{T}[\bm{w}\odot(\bm{X}\bm{u})]$. Therefore, we have

Furthermore, we know that the Hessian matrices of $f(\bm{u},\bm{v},\bm{w})$ with respect to $\bm{u}$, $\bm{v}$ and $\bm{w}$ are zero matrix, respectively. Thus, $\nabla_{\bm{u}}f(\bm{u},\bm{v},\bm{w})$, $\nabla_{\bm{v}}f(\bm{u},\bm{v},\bm{w})$ and $\nabla_{\bm{w}}f(\bm{u},\bm{v},\bm{w})$ are Lipschitz continuous, and their Lipschitz constants can be any constant greater than zero ($L_{u}>0$, $L_{v}>0$, and $L_{w}>0$).

Recently, the PALM algorithm [18] has been proposed to solve a class of non-convex and non-smooth problems. Based on the PALM framework, we develop a block proximal gradient algorithm to solve (4) by iteratively updating each variable when fixing others.

1) Computation of $\bm{u}^{t+1}$ using $\bm{u}^{t}$, $\bm{v}^{t}$ and $\bm{w}^{t}$. We perform the proximal gradient step for $\bm{u}^{t+1}$ when $\bm{v}$ and $\bm{w}$ are fixed:

where $\overline{\bm{u}}=\bm{u}^{t}-\frac{1}{L_{u}}\nabla_{\bm{u}}f(\bm{u}^{t},\bm{v}^{t},\bm{w}^{t})$. To solve the above problem (6), we define an $\ell_{0}$-ball projection problem as Eq. (7) and give its closed-form solution in Proposition 1.

Based on Proposition 1, we obtain the closed-form solution of Eq. (6) as follows:

where $\overline{\bm{u}}=\bm{u}^{t}-\frac{1}{L_{u}}\nabla_{\bm{u}}f(\bm{u}^{t},\bm{v}^{t},\bm{w}^{t})$, the partial derivatives $\nabla_{\bm{u}}f(\bm{u},\bm{v},\bm{w})=-\bm{X}^{T}[\bm{w}\odot(\bm{Y}\bm{v})]$ and the Lipschitz constant $L_{u}$ can be set to any constant greater than zero.

2) Computation of $\bm{v}^{t+1}$ using $\bm{u}^{t+1}$, $\bm{v}^{t}$ and $\bm{w}^{t}$. We perform the proximal gradient step for $\bm{v}^{t+1}$ when $\bm{u}$ and $\bm{w}$ are fixed:

where $\overline{\bm{v}}=\bm{v}^{t}-\frac{1}{L_{v}}\nabla_{\bm{u}}f(\bm{u}^{t+1},\bm{v}^{t},\bm{w}^{t})$. Based on the Proposition 1, we obtain the closed-form solution of Eq. (9) as follows:

where $\overline{\bm{v}}=\bm{v}^{t}-\frac{1}{L_{v}}\nabla_{\bm{u}}f(\bm{u}^{t+1},\bm{v}^{t},\bm{w}^{t})$, the partial derivatives $\nabla_{\bm{v}}f(\bm{u},\bm{v},\bm{w})=-\bm{Y}^{T}[\bm{w}\odot(\bm{X}\bm{u})]$ and the Lipschitz constant $L_{v}$ can be set to any constant greater than zero.

3) Computation of $\bm{w}^{t+1}$ using $\bm{u}^{t+1}$, $\bm{v}^{t+1}$ and $\bm{w}^{t}$. We perform the proximal gradient step for $\bm{w}^{t+1}$ when $\bm{u}$ and $\bm{v}$ are fixed:

where $\overline{\bm{w}}=\bm{w}^{t}-\frac{1}{L_{w}}\nabla_{\bm{w}}f(\bm{u}^{t+1},\bm{v}^{t+1},\bm{w}^{t})$. To solve the above model (11), we define a projection function $\widehat{\bm{w}}:=P_{\infty,+}(\bm{w})$ for an arbitrary vector $\bm{w}$, which satisfies as follows for all $j$:

We obtain the closed-form solution of Eq. (11) using the following Proposition 2.

Based on Proposition 2, we obtain the update rule of $\bm{w}^{t+1}$:

where $\overline{\bm{w}}=\bm{w}^{t}-\frac{1}{L_{w}}\nabla_{\bm{w}}f(\bm{u}^{t+1},\bm{v}^{t+1},\bm{w}^{t})$, $P_{\infty,+}(\cdot)$ is defined in Eq. (12), $\nabla_{\bm{w}}f(\bm{u},\bm{v},\bm{w})=-(\bm{X}\bm{u})\odot(\bm{Y}\bm{v})$ and the Lipschitz constant $L_{w}$ can be set to any constant greater than zero.

Combing (8), (10) and (13), we propose the block proximal gradient algorithm to solve Eq. (4) (See Algorithm 1).

Initialization. We assume that all the samples have equal contribution and initialize $\bm{w}$ as $\bm{w}^{0}=\bm{1}$ (with all entries are ones) in Algorithm 1. Furthermore, Algorithm 1 is run with multiple random initial points $\bm{u}^{0}$ and $\bm{v}^{0}$. Finally, the model with the optimal objective function value is selected.

Stopping Condition. We can use the following ways to stop the algorithm: (i) the stopping criterion can be set as follows:

where $\bm{u}^{t}$, $\bm{v}^{t}$, and $\bm{w}^{t}$ represent the $t$-th iteration of $\bm{u}$, $\bm{v}$, and $\bm{w}$ respectively. (ii) the stopping criteria can be set as follows:

where $\mbox{obj}(t)$ is the objective function value after the $t$-th iteration.

The computational complexity of Algorithm 1 mainly depends on the computation of Eq. (8), Eq. (10) and Eq. (13). Computing Eq. (8) is mainly composed of the computation of $\nabla_{\bm{u}}f(\bm{u},\bm{v},\bm{w})=-\bm{X}^{T}[\bm{w}\odot(\bm{Y}\bm{v})]$ which consumes $\mathcal{O}(np+nq+n)$. Computing Eq. (10) is mainly composed of the computation of $\nabla_{\bm{v}}f(\bm{u},\bm{v},\bm{w})=-\bm{Y}^{T}[\bm{w}\odot(\bm{X}\bm{u})]$ which consumes $\mathcal{O}(nq+np+n)$. Similarly, Computing Eq. (13) is mainly composed of the computation of $\nabla_{\bm{w}}f(\bm{u},\bm{v},\bm{w})=-(\bm{X}\bm{u})\odot(\bm{Y}\bm{v})$ which consumes $\mathcal{O}(nq+np+n)$. Overall, the complexity of Algorithm 1 is $\mathcal{O}(tnp+tnq+tn)$, where $t$ is the number of iterations and it is empirically set to 20 in our experiments.

In this subsection, we show that Algorithm 1 has theoretical convergence guarantees (Theorem 1).

A co-module is decide with two steps: (1) extract a submatrix in the data marix $\bm{X}\in\mathbb{R}^{n\times p}$, where rows and columns are corresponding to non-zero entries in $\bm{w}$ and $\bm{u}$, respectively; (2) extract another submatrix in the data matrix $\bm{Y}\in\mathbb{Y}^{n\times q}$, where rows and columns are corresponding to non-zero entries in $\bm{w}$ and $\bm{v}$, respectively. The combination of two extracted sub-matrices is called a co-module (see Figure 1).

We explain how wsPLS identify multiple co-modules. After identifying the first co-module using Algorithm 1, we update matrices $\bm{X}$ and $\bm{Y}$ by removing their rows (samples) in the first co-module, and we repeatedly apply Algorithm 1 to identify the next co-module on the updated data matrices (see the three co-modules in the rightmost subgraph of Figure 1).

Assume that multiple matrices $\bm{X}_{i}\in\mathbb{R}^{n\times p_{i}}$ ($i=1,\cdots,m$) across a common sample set are given and their columns are normalized with zero-mean and unit-variance. Similar to [32], we propose the multi-view wsPLS (mwsPLS) model using a sum way, which is called mwsPLS-scheme1 for short as follows:

Eq. (16) is equivalent to

Optimization. We observe that the objective function of Eq. (17) satisfies $f(\bm{u},\bm{w})=-\sum\limits_{i<j}\bm{w}^{T}[(\bm{X}_{i}\bm{u}_{i})\odot(\bm{X}_{j}\bm{u}_{j})]=-\bm{u}_{i}^{T}\bm{X}_{i}^{T}[\sum\limits_{i\neq j}\bm{w}\odot(\bm{X}_{j}\bm{u}_{j})]+C$, where $C$ is independent for any $\bm{u}_{i}$. So, we can calculate the partial derivatives of $\bm{u}_{i}$ ($i=1,\cdots,m$) and $\bm{w}$ as follows:

Similar to Algorithm 1, we develop the block proximal gradient algorithm to solve (17) by using iteratively updating way. The detail algorithm procedure is shown in Algorithm 2 where Eq. (19) and Eq. (20) are solved using Proposition 1 and Proposition 2, respectively.

Initialization. We set $\bm{w}^{0}=\bm{1}$ (all entries are ones) in Algorithm 2 and randomly initialize {$\bm{u}_{i}^{0},\leavevmode\nobreak\ i=1,\cdots,m$} from a normal distribution. Algorithm 2 is run with multiple random initial points and the model with the optimal value is selected.

Complexity Analysis. The computational burden of Algorithm 2 mainly depends on the computation of Eq. (19) and Eq. (20). The computation of Eq. (19) is mainly composed of the computation of $\nabla_{\bm{u}_{i}}f(\bm{u},\bm{w})=-\bm{X}_{i}^{T}\sum_{j\neq i}[\bm{w}\odot(\bm{X}_{j}\bm{u}_{j})]$ which consumes $\mathcal{O}(\sum_{i}np_{i})$. The computation of Eq. (20) is mainly composed of the computation of $\nabla_{\bm{w}}f(\bm{u},\bm{w})=-\sum_{i<j}[(\bm{X}_{i}\bm{u}_{i})\odot(\bm{X}_{j}\bm{u}_{j})]$ which consumes $\mathcal{O}((m-1)\sum_{i}np_{i})$ because $\sum_{i<j}(np_{i}+np_{j})=(m-1)\sum_{i}np_{i}$. In total, the complexity of Algorithm 2 is $\mathcal{O}(tmnp_{1}+tmnp_{2}+\cdots+tmnp_{m})$, where $t$ is the number of iterations and it is empirically set to 20 in our experiments.

Convergence Analysis. Similar to Lemma 1, we can prove that the objective function of the mwsPLS-scheme1 model satisfies the KŁ property. Therefore, Algorithm 2 also has theoretical convergence.

We observe $\bm{u}_{1}^{T}\bm{X}_{1}^{T}\mbox{diag}(\bm{w})\bm{X}_{2}\bm{u}_{2}=\bm{w}^{T}[(\bm{X}_{1}\bm{u}_{1})\odot(\bm{X}_{2}\bm{u}_{2})]$, so that we intuitively extend wsPLS for multi-view data analysis using a product way. Herein, we propose the second mwsPLS model using the product way, which is called mwsPLS-scheme2 for short as follows:

where $\bigodot\limits_{i=1}^{m}(\bm{X}_{i}\bm{u}_{i})=(\bm{X}_{1}\bm{u}_{1})\odot(\bm{X}_{2}\bm{u}_{2})\odot\cdots\odot(\bm{X}_{m}\bm{u}_{m})$.

Optimization. Formally, the objective function in Eq. (21) is equivalent to

where $\bm{z}_{i}=\bigodot\limits_{j\neq i}^{m}(\bm{X}_{j}\bm{u}_{j})$. Therefore, we calculate the partial derivatives for $\bm{u}_{i}$ ($i=1,\cdots,m$) and $\bm{w}$ as follows:

Finally, we solve the optimization problem (21) using a similar manner as Algorithm 2 and the detail procedure is shown in Algorithm 3.

Initialization. The initial point setting of Algorithm 3 is the same as that of Algorithm 2, and the description is omitted here.

Complexity Analysis. The complexity of $\nabla_{\bm{u}_{i}}f(\bm{u},\bm{w})=-\bm{X}_{i}^{T}(\bm{w}\odot\bm{z}_{i})$ and $\bm{z}_{i}=\bigodot\limits_{j\neq i}^{m}(\bm{X}_{j}\bm{u}_{j})$ is $\mathcal{O}(\sum_{i}np_{i})$. The complexity of $\nabla_{\bm{w}}f(\bm{u},\bm{w})=-\bigodot\limits_{i=1}^{m}(\bm{X}_{i}\bm{u}_{i})$ is also $\mathcal{O}(\sum_{i}np_{i})$. So, the complexity of Algorithm 3 is $\mathcal{O}(tmnp_{1}+tmnp_{2}+\cdots+tmnp_{m})$, where $t$ is the number of iterations and it is empirically set to 20 in our experiments.

Convergence Analysis. Similar to Lemma 1, we can also prove that the objective function of the mwsPLS-scheme2 model satisfies the KŁ property. So, Algorithm 3 also has theoretical convergence guarantees and globally converges to a critical point.

We compare the performance of the proposed methods with PLS and its variants:

• •

  PLS: It is a classical model without any sparse constrain [29]:

  		$\displaystyle\underset{\bm{u},\bm{v}}{\text{maximize}}$		$\displaystyle\bm{u}^{T}\bm{X}^{T}\bm{Y}\bm{v}$
  		subject to		$\displaystyle\|\bm{u}\|=\|\bm{v}\|=1,$

  which is solved using the single value decomposition (SVD) method for $\bm{X}^{T}\bm{Y}$.

• •

  PMD-sPLS: It is a sparse model which uses the $\ell_{1}$-norm constraint, which is inspired by [32]:

  		$\displaystyle\underset{\bm{u},\bm{v}}{\text{maximize}}$		$\displaystyle\bm{u}^{T}\bm{X}^{T}\bm{Y}\bm{v}$
  		subject to		$\displaystyle\|\bm{u}\|=\|\bm{v}\|=1,\|\bm{u}\|<c_{1},\|\bm{v}\|<c_{2},$

  which is solved by the penalized matrix decomposition (PMD) method [32].

• •

  $\ell_{0}$-sPLS: It is a sparse PLS model which uses the $\ell_{0}$-norm constraint, which is inspired by [30, 33]:

  		$\displaystyle\underset{\bm{u},\bm{v}}{\text{maximize}}$		$\displaystyle\bm{u}^{T}\bm{X}^{T}\bm{Y}\bm{v}$
  		subject to		$\displaystyle\|\bm{u}\|=\|\bm{v}\|=1,\|\bm{u}\|_{0}\leq k_{u},\|\bm{v}\|_{0}\leq k_{v},$

  which is solved using an alternating iterative algorithm.

• •

  $\ell_{2}/\ell_{0}$-wsPLS: It is a weighted sparse PLS model which uses the $\ell_{2}/\ell_{0}$-norm constrains for the weighed vector, which is inspired by [31].

  		$\displaystyle\underset{\bm{u},\bm{v},\bm{w}}{\text{maximize}}$		$\displaystyle\bm{u}^{T}\bm{X}^{T}\mbox{diag}(\bm{w})\bm{Y}\bm{v}$
  		subject to		$\displaystyle\|\bm{u}\|_{0}\leq k_{u},\|\bm{v}\|_{0}\leq k_{v},\|\bm{u}\|=\|\bm{v}\|=1,$
  		$\displaystyle\|\bm{w}\|_{2}=1,\|\bm{w}\|_{0}\leq k_{w},$

  which is also solved using the proposed Algorithm 1 with a small change.

To generate the synthetic data of $\bm{X}\in\mathbb{R}^{n\times p}$ and $\bm{Y}\in\mathbb{R}^{n\times q}$ ($n$ is the number of samples, $p$ and $q$ are the number of two different features, respectively), we firstly generated $\bm{u}$, $\bm{v}$ and $\bm{w}$, and then generate them by

where the elements of $\bm{\epsilon}_{1}$ and $\bm{\epsilon}_{2}$ are randomly sampled from a standard normal distribution and $\gamma_{1}$, $\gamma_{2}$ are two nonnegative parameters to control the signal-to-noise ratio (SNR) which is defined by

where $\mathbf{E}(\|\gamma_{1}\bm{\epsilon}_{1}\|_{F}^{2})$ and $\mathbf{E}(\|\gamma_{2}\bm{\epsilon}_{2}\|_{F}^{2})$ denote the expected sum of squares of noise.

(1) Simulation data I ($n$ = 50, $p$ = 80, $p$ = 100). We firstly generate three vectors $\bm{u}$, $\bm{v}$ and $\bm{w}$ using

where $rep(x,n)$ denotes a row vector of size $n$ with all elements are equal to $x$. We then generate the simulation data I using the formula (26) with $\mbox{SNR}_{1}=\mbox{SNR}_{2}=0.1$.

(2) Simulation data II ($n$ = 100, $p$ = 800, $p$ = 1000). We firstly generate three vectors $\bm{u}$, $\bm{v}$ and $\bm{w}$ using