Simultaneous sparse estimation of canonical vectors in the p>> N setting

Paper can be found in
http://arxiv.org/abs/1403.6095

1. Goal

In this paper, the authors are concerned with multi-group classification(LDA) and want to estimate G-1 canonical vectors simultaneously rather than estimate in a sequential fashion.

2. Method

2.1 ratio-trace version of LDA

Suppose $\Sigma_b \in \mathbb{R}^{p\times p}$ and $\Sigma_w\in\mathbb{R}^{p\times p}$ are population between and within group covariance matrices. Note that ratio-trace version of linear discriminant analysis is dealing with the following optimization problem:

$$\begin{align} \max\quad&\text{Tr}\left(\left(V^T\Sigma_w V\right)^{-1}V^T\Sigma_b V\right)\\ \end{align}$$
which can be solved via generalized eigen decomposition of $\Sigma_w^{-1}\Sigma_b$.

proof: Taking derivative with respect to $V$ and setting to zero,

$$\begin{align} &\frac{\partial \text{Tr}\left(\left(V^T\Sigma_w V\right)^{-1}V^T\Sigma_b V\right)}{\partial V}\\ &=-2\Sigma_w V \left(V^T\Sigma_w V\right)^{-1}V^T\Sigma_b V\left(V^T\Sigma_w V\right)^{-1}+2\Sigma_b V\left(V^T\Sigma_w V\right)^{-1}=0 \end{align}$$

how to get this formula can be found in
http://blog.youkuaiyun.com/comeyan/article/details/50514610

we have

$$\begin{align} &\Sigma_w V \left(V^T\Sigma_w V\right)^{-1}V^T\Sigma_b V = \Sigma_b V\\ \Rightarrow &\Sigma_w^{-1} \Sigma_b V=V \left(V^T\Sigma_w V\right)^{-1}V^T\Sigma_b V \\ \Rightarrow&\Sigma_w^{-1} \Sigma_b V = V B \Theta B^{-1}\\ \Rightarrow&\Sigma_w^{-1} \Sigma_b VB = VB\Theta \end{align}$$

Note that:

• The third equivalence is because we can simultaneously diagonalize $V^T\Sigma_w V$ and $V^T\Sigma_ V$ using $B$, which is to say($\Theta$ is a diagonal matrix)
  $$B^{-1}V^T\Sigma_w V B = I, B^{-1}V^T\Sigma_b VB = \Theta$$
  So
  $$V^T\Sigma_w V = I, V^T\Sigma_b V = B\Theta B^{-1}$$
  how to diagonalize two matrices simultaneously can be found in
  http://blog.youkuaiyun.com/comeyan/article/details/50521034
• The fourth equivalence tells us that $VB$ is the generalized eigen vectors of $\Sigma_w^{-1}\Sigma_b$

Summarize:
From the above analysis, we know that to solve the ratio-trace version of LDA is to solve a generalized eigen decomposition problem, which can be done by diagonaling these two matrices simultaneously.

2.2 Model

Since the ratio-trace version of LDA is to solve a generalized eigen decomposition problem of $\Sigma_w^{-1}\Sigma_b$ and the eigenvectors are unique only up to normalization, Taking advantage of the unique of the eigenspace in defining a scale-invariant classification rule.

Notations to be used:

• $G$ is the total number of groups.
• $\text{rank}(\Sigma_b)=G-1$
• $\mu_i, i=1,2,\cdots, G$ is the mean of each group.
• $\pi_i, i=1,2,\cdots,G$ is the prior probability of each group

The goal of this paper is to find find $G-1$ eigenvectors $\Phi$ corresponding to non-zero eigenvalues of $\Sigma_w^{-1}\Sigma_b$. *

Note that the $G-1$ vectors can be formularized in a closed form.
Proposition 1(population version). The following decomposition holds: $\Sigma_b = \Delta\Delta^T$, where for $r=1,2,\cdots, G-1$ the $r$th column of $\Delta$ has the form

$$\Delta_r = \frac{\sqrt{\pi_{r+1}}\left(\sum_{i=1}^r\pi_i(\mu_i - \mu_{r+1})\right)}{\sqrt{\sum_{i=1}^r \pi_r \sum_{i=1}^{r+1}\pi_i}}$$

Proposition 2(sample version). The following decomposition holds: $\hat\Sigma_b =DD^T$, where for $r=1,2,\cdots, G-1$ the $r$th column of $D$ has the form

$$D_r = \frac{\sqrt{n_{r+1}}\left(\sum_{i=1}^rn_i(\hat \mu_i - \hat\mu_{r+1})\right)}{\sqrt{N}\sqrt{\sum_{i=1}^r n_i \sum_{i=1}^{r+1}n_i}}$$

the proof of sample version can be found in(using orthogonal contrasts of unbalanced data)
http://blog.youkuaiyun.com/COMEYAN/article/details/50521276
orthogonal contrasts can be found in

Formularize $\Sigma_b$ in this form of low-rank decomposition is because it has a closed form and have intuitive interpretation in terms of the differences between the group means.

Then give the closed form of generalized eigen vectors $\Phi$

Proposition 3. Define $\Delta$ and $D$ as in last two propositions. There exists a matrix $P\in \mathbb{O}^{G-1}$ such that

$$\Phi = \Sigma_w^{-1}\Delta P$$
Moreover, if $\hat\Sigma_w$ is nonsingular, there exists a matrix $R\in \mathbb{O}^{G-1}$ such that
$$\hat\Phi = \hat\Sigma_w^{-1}D R$$

If we use Mahalanobis distance for classification, the classification function is free of orthogonal scale, so we can use a simple projection matrix instead.

$$\begin{align} \tilde \Phi &= \Sigma_w^{-1}\Delta\\ \end{align}$$

Now it is sufficient to estimate $\tilde \Phi = \Sigma_w^{-1}\Delta$, which can be defined as

$$\tilde\Phi = \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\left\|\Sigma_w^{1/2} V - \Sigma_w^{-1/2}\Delta\right\|_F^2 = \frac{1}{2}\text{tr}\left(V^T\Sigma_w V - 2V^T \Delta\right)$$

Sample-version
Model step 1

$$\tilde V = \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\hat\Sigma_w V - 2V^TD\right)$$

To get sparse solution, adding group-lasso penalty
Model step 2

$$\hat V = \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\hat\Sigma_w V - 2V^TD\right) +\lambda\sum_{i=2}^p \|V_{i*}\|_2$$

But the objective function can be unbounded below when $\hat\Sigma_w$ is singular. So regularization of $\hat\Sigma_w$ is needed, $\tilde \Sigma_w = \hat\Sigma_w+\rho I$

Model step 3

$$\begin{align} \hat V &= \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\tilde\Sigma_w V - 2V^TD\right) +\lambda\sum_{i=2}^p \|V_{i*}\|_2 \\ &= \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\hat\Sigma_w V \right) + \frac{\rho}{2}\left\|V - \frac{1}{\rho}D\right\|_F^2 + \lambda\sum_{i=2}^p \|V_{i*}\|_2 \end{align}$$

But this means letting $\hat V$ being close to $\frac{1}{\rho}D$, which is not true. So replacing $\hat\Sigma_w$ by $\hat\Sigma_t$, we have the final model

Model step 4

$$\begin{align} \hat V &= \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\hat\Sigma_t V - 2V^TD\right) +\lambda\sum_{i=2}^p \|V_{i*}\|_2 \\ &= \arg\min_{V\in\mathbb{R}^{p\times G-1}} \frac{1}{2}\text{tr}\left(V^T\hat\Sigma_w V \right) + \frac{1}{2}\left\|V^TD - I\right\|_F^2 + \lambda\sum_{i=2}^p \|V_{i*}\|_2 \end{align}$$

• the first item is to minimize the within group variability,
• the second item is to control the level of the between group variability
• the third item is to induce sparsity

3. Theory

The model has model selection property and the misclassification error coincides with the population rule.

4. Algorithm

Block coordinate descent algorithm is used to get the solution. As it takes advantage of warm starts when solving for a range of tuning parameters and is onr of the fastest algorithm for smooth losses with separable regularizers.