 Research
 Open Access
 Published:
Sparse representation approaches for the classification of highdimensional biological data
BMC Systems Biologyvolume 7, Article number: S6 (2013)
Abstract
Background
Highthroughput genomic and proteomic data have important applications in medicine including prevention, diagnosis, treatment, and prognosis of diseases, and molecular biology, for example pathway identification. Many of such applications can be formulated to classification and dimension reduction problems in machine learning. There are computationally challenging issues with regards to accurately classifying such data, and which due to dimensionality, noise and redundancy, to name a few. The principle of sparse representation has been applied to analyzing highdimensional biological data within the frameworks of clustering, classification, and dimension reduction approaches. However, the existing sparse representation methods are inefficient. The kernel extensions are not well addressed either. Moreover, the sparse representation techniques have not been comprehensively studied yet in bioinformatics.
Results
In this paper, a Bayesian treatment is presented on sparse representations. Various sparse coding and dictionary learning models are discussed. We propose fast parallel activeset optimization algorithm for each model. Kernel versions are devised based on their dimensionfree property. These models are applied for classifying highdimensional biological data.
Conclusions
In our experiment, we compared our models with other methods on both accuracy and computing time. It is shown that our models can achieve satisfactory accuracy, and their performance are very efficient.
Introduction
The studies in biology and medicine have been revolutionarily changed since the invents of many highthroughput sensory techniques. Using these techniques, the molecular phenomenons can be probed with a high resolution. In the virtue of such techniques, we are able to conduct systematic genomewide analysis. In the last decade, many important results have been achieved by analyzing the highthroughput data, such as microarray gene expression profiles, gene copy numbers profiles, proteomic mass spectrometry data, nextgeneration sequences, and so on.
On one hand, biologists are enjoining the richness of their data; one another hand, bioinformaticians are being challenged by the issues of the highdimensional data. Many of the analysis can be formulated into machine learning tasks. First of all, we have to face to the cures of high dimensionality which means that many machine learning models can be overfitted and therefore have poor capability of generalization. Second, if the learning of a model is sensitive to the dimensionality, the learning procedure could be extremely slow. Third, many of the data are very noise, therefore the robustness of a model is necessary. Forth, the highthroughput data exhibit a large variability and redundancy, which make the mining of useful knowledge difficult. Moreover, the observed data usually do not tell us the key points of the story. We need to discover and interpret the latent factors which drive the observed data.
Many of such analysis are classification problem from the machine learning viewpoint. Therefore in this paper, we focus our study on the classification techniques for highdimensional biological data. The machine learning techniques addressing the challenges above can be categorized into two classes. The first one aims to directly classify the highdimensional data while keeping the good capability of generalization and efficiency in optimization. The most popular method in this class is the regularized basisexpended linear model. One example is the stateoftheart support vector machine (SVM) [1]. SVM can be kernelized and its result is theoretically sound. Combining different regularization terms and various loss functions, we can have many variants of such linear models [2]. In addition to classification, some of the models can be applied to regression and feature (biomarker) identification. However, most of the learned linear models are not interpretable, while interpretability is usually the requirement of biological data analysis. Moreover, linear models can not be extended naturally to multiclass data, while in bioinformatics a class may be composed of many subtypes.
Another technique of tackling with the challenges is dimension reduction which includes feature extraction and feature selection. Principal component analysis (PCA) [3] is the oldest feature extraction method and is widely used in processing highdimensional biological data. The basis vectors produced by PCA are orthogonal, however many patterns in bioinformatics are not orthogonal at all. The classic factor analysis (FA) [4] also has orthogonal constraints on the basis vectors, however its Bayesian treatment does not necessarily produce orthogonal basis vectors. Bayesian factor analysis will be introduced in the next section.
Sparse representation (SR) [5] is a parsimonious principle that a sample can be approximated by a sparse linear combination of basis vectors. Nonorthogonal basis vectors can be learned by SR, and the basis vectors may be allowed to be redundant. SR highlights the parsimony in representation learning [6]. This simple principle has many strengthes that encourage us to explore its usefulness in bioinformatics. First, it is very robust to redundancy, because it only select few among all of the basis vectors. Second, it is very robust to noise [7]. Furthermore, its basis vectors are nonorthogonal, and sometimes are interpretable due to its sparseness [8]. There are two techniques in SR. First, given a basis matrix, learning the sparse coefficient of a new sample is called sparse coding. Second, given training data, learning the basis vector is called dictionary learning. As dictionary learning is, in essence, a sparse matrix factorization technique, nonnegative matrix factorization (NMF) [9] can be viewed a specific case of SR. For understanding sparse representation better, we will give the formal mathematical formulation from a Bayesian perspective in the next section.
This paper is the significant extension of our preliminary work presented in [10] where sparse representation is treated from regularization and optimization perspectives. In this paper, we formulate sparse representation from a Bayesian viewpoint. We show that using different prior distributions, we can obtain various sparse coding and dictionary learning models. Although there exists some works, for example [11], which apply sparse coding in the classification of biological data, to the best of our knowledge, this is the first time that sparse representation is intensively and systematically studied in the area of bioinformatics. This study has the following contributions:

1.
We give a Bayesian treatment on the sparse representation, which is very helpful to understand and design sparse representation models.

2.
Kernel sparse coding techniques are proposed for direct classification of highdimensional biological data.

3.
Fast parallel activeset algorithms are devised for sparse coding.

4.
An efficient generic framework of kernel dictionary learning for feature extraction is proposed.

5.
We reveal that the optimization and decision making in sparse representation is dimensionfree.
We organize the rest of this paper as follow. We first introduce factor analysis and sparse representation from a Bayesian aspect. Classification method based on sparse coding is then introduced and the activeset methods are proposed for the optimization. Their kernel extensions are given as well. Then various dictionary learning models are given. After that, a generic optimization framework is devised to optimize these models. In the same section, dictionarylearningbased classification and its kernel extension are proposed as well. Then we describe our computational experiments on two highdimensional data sets. Finally, conclusions and future works are drawn.
Related work from a Bayesian viewpoint
Both (sparse) factor analysis and sparse representation models can be used as dimension reduction techniques. Due to their intuitive similarity, it is necessary to give their definitions for comparison. In this section, we briefly survey the sparse factor analysis and sparse representation in a Bayesian viewpoint. The introduction of sparse representation is helpful to understand the content of the subsequent sections. Hereafter, we use the following notations unless otherwise stated. Suppose the training data is $D\in {\mathbb{R}}^{m\times n}$ (m is the number of features and n is the number of training instances (samples or observations)), the class information is in the ndimensional vector c. Suppose p new instances are in $B\in {\mathbb{R}}^{m\phantom{\rule{2.77695pt}{0ex}}\times p}$.
Sparse Bayesian (latent) factor analysis
The advantages of Bayesian (latent) factor analysis model [12] over likelihood (latent) factor analysis model are that

1.
The knowledge regarding the model parameters from experts and previous investigations can be incorporated through the prior.

2.
The values of parameters are refined using the current training observations.
The Bayesian factor analysis model [12] can be formulated as
where $b\in {\mathbb{R}}^{m\phantom{\rule{2.77695pt}{0ex}}\times \phantom{\rule{0.3em}{0ex}}1}$ is an observed multivariate variable, $\mu \in {\mathbb{R}}^{m\phantom{\rule{2.77695pt}{0ex}}\times \phantom{\rule{0.3em}{0ex}}1}$ is the population mean, $A\in {\mathbb{R}}^{m\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\times \phantom{\rule{0.3em}{0ex}}k}$ is latent factor loading matrix, and $x\in {\mathbb{R}}^{k\phantom{\rule{0.3em}{0ex}}\times \phantom{\rule{0.3em}{0ex}}1}$ is latent factor score (k ≪ m), and $\epsilon \in {\mathbb{R}}^{m\phantom{\rule{0.3em}{0ex}}\times \phantom{\rule{0.3em}{0ex}}1}$ is an idiosyncratic error term. This model is restricted by the following constraints or assumptions:

1.
The error term is normally distributed with mean 0 and covariance Φ: ε~ N (0, Φ). Φ is diagonal on average.

2.
The factor score vector is also normally distributed with mean 0 and identity covariance R= I: x~ N (0, R); and the factor loading vector is normally distributed: a_{ i }~ N (0, Δ) where Δ is diagonal. Alternatively, the factor loading vectors can be normally distributed with mean 0 and identity covariance Δ = I; and the factor score vector is normally distributed with mean 0 and diagonal covariance R. The benefit of identity covariance either on x or A is that arbitrary scale interchange between A and x due to scale invariance can be avoided.

3.
x is independent of ε.
For n training instances D, we have the likelihood:
where tr(M) is the trace of square matrix M.
The variants of Bayesian factor analysis models differ in the decomposition of the joint priors. The simplest one may be p(µ, A , Y) = p(µ)p(A)p(Y). Suppose k is fixed a priori. The posterior therefore becomes
The model parameters are usually estimated via Markov chain Monte Carlo (MCMC) techniques.
Sparse Bayesian factor analysis model imposes a sparsityinducing distribution over the factor loading matrix instead of Gaussian distribution. In [13], the following mixture of prior is proposed:
where π_{ ij } is the probability of a nonzero a_{ ij } and δ_{0}(·) is the Dirac delta function at 0. Meanwhile, A is constrained using the lower triangular method. Bayesian factor regression model (BFRM) is the combination of Bayesian factor analysis and Bayesian regression [13]. It has been applied in oncogenic pathway studies [4] as a variable selection method.
Sparse representation
Sparse representation (SR) is a principle that a signal can be approximated by a sparse linear combination of dictionary atoms [14]. The SR model can be formulated as
where $A=\left[{a}_{1},\cdots \phantom{\rule{0.3em}{0ex}},{a}_{k}\right]$ is called dictionary, a_{ i }is a dictionary atom, x is a sparse coefficient vector, and ε is an error term. A, x, and k are the model parameters. SR model has the following constraints:

1.
The error term is Gaussian distributed with mean zero and isotropic covariance, that is ε~ N (0, Φ) where Φ = ϕ I where ϕ is a positive scalar.

2.
The dictionary atoms is usually Gaussian distributed, that is a_{i} ~ N (0, Δ) where Δ = I. The coefficient vector should follows a sparsityinducing distribution.

3.
x is independent of ε.
Through comparing the concepts of Bayesian factor analysis and Bayesian sparse representation, we can find that the main difference between them is that the former applies a sparsityinducing distribution over the factor loading matrix, while the later uses a sparsityinducing distribution on the factor score vector.
Sparse representation involves sparse coding and dictionary learning. Given a new signal b and a dictionary A, learning the sparse coefficient x is termed sparse coding. It can be statistically formulated as
Suppose the coefficient vector has Laplacian prior with zero mean and isotropic variance, that is $p\left(x\mathrm{\Gamma}\right)=L\left(0,\mathrm{\Gamma}\right)=\frac{1}{{\left(2\gamma \right)}^{k}}{e}^{\frac{{x}_{1}}{\gamma}}$. The likelihood is Gaussian distributed as $p\left(bA,x,\mathrm{\Phi}\right)=N\left(Ax,\mathrm{\Phi}\right)=\frac{1}{{\left(2\pi \right)}^{\frac{m}{2}}{\varphi}^{\frac{m}{2}}}{e}^{\frac{1}{2\varphi}{bAx}_{2}^{2}}$. The posterior is thus
Taking the logarithm, the above is thus
where c is a constant term. We can see that maximizing the posterior is equivalent to minimizing the following task:
where $\lambda =\frac{\varphi}{\gamma}$. Hereafter we call Equation (10) l_{1}least squares (l_{1} LS) sparse coding model. It is known as the l_{1}regularized regression model in regularization theory. It coincides with the wellknown LASSO model [15], which in fact is a maximum a posteriori (MAP) estimation.
Given training data D, learning (or estimating) the dictionary A, the coefficient vectors Y, and the number of dictionary atoms k is called dictionary learning. Suppose k is given a priori, and consider the Laplacian prior over Y and the Gaussian prior over A, and suppose $p\left(A,Y\right)=p\left(A\right)p\left(Y\right)=\prod _{i=1}^{k}\left(p\left({a}_{i}\right)\right)\prod _{i=1}^{n}\left(p\left({y}_{i}\right)\right)$. We thus have the prior:
The likelihood is
The posterior is
Ignoring the normalization term (that is the marginal likelihood), the logposterior is thus
Therefore the MAP estimation of dictionary learning task is
where α = f and $\lambda =\frac{\varphi}{\gamma}$. Equation (16) is known as a dictionary learning model based on l_{1}regularized least squares.
In the literature, the kernel extension of sparse representation is realized in two ways. The first way is to use empirical kernel in sparse coding as in [16], where dictionary learning is not considered. The second way is the one proposed in [17], where dictionary learning is involved. However, the dictionary atoms are represented and updated explicitly. This could be intractable, as the number of dimensions of dictionary atoms in the feature space is very high even infinite. In the later sections, we give our kernel extensions of sparse coding and dictionary learning, respectively, where any kernel functions can be used and the dictionary is updated efficiently.
Sparse coding methods
Since, the l_{1}LS sparse coding (Equation (10)) is a twosided symmetric model, thus a coefficient can be zero, positive, or negative [18]. In Bioinformatics, l_{1}LS sparse coding has been applied for the classification of microarray gene expression data in [11]. The main idea is in the following. First, training instances are collected in a dictionary. Then, a new instance is regressed by l_{1}LS sparse coding. Thus its corresponding sparse coefficient vector is obtained. Next, the regression residual of this instance to each class is computed, and finally this instance is assigned to the class with the minimum residual.
We generalize this methodology in the way that the sparse code can be obtained by many other regularization and constraints. For example, we can pool all training instances in a dictionary (hence k = n and A= D), and then learn the nonnegative coefficient vectors of a new instance, which is formulated as an onesided model:
We called this model the nonnegative least squares (NNLS) sparse coding. NNLS has two advantages over l_{1}LS. First, the nonnegative coefficient vector is more easily interpretable than coefficient vector of mixed signs, under some circumstances. Second, NNLS is a nonparametric model. From a Bayesian viewpoint, Equation (17) is equivalent to the MAP estimation with the same Gaussian error as in Equation (6), but the following discrete prior:
This nonnegative prior implies that, the elements in x are independent, and the probability of x_{ i } = 0 is 0.5 and the probability of x_{ i } > 0 is 0.5 as well. (That is the probabilities of x_{ i } being either 0 or positive are equal, and the probability of being negative is zero.) Inspired by many sparse NMFs, l_{1}regularization can be additionally used to produce more sparse coefficients than NNLS above. The combination of l_{1}regularization and nonnegativity results in the l_{1}NNLS sparse coding model as formulated below:
We call Equation (19) the l_{1}NNLS model. It is more flexible than NNLS, because it can produce more sparse coefficients as controlled by λ. This model in fact uses the following prior:
Now, we give the generalized sparsecodingbased classification approach in details. The method is depicted in Algorithm 1. We shall give the optimization algorithms, later, required in the first step. The NN rule mentioned in Algorithm 1 is inspired by the usual way of using NMF as a clustering method. Suppose there are C classes with labels 1, ⋯, C. For a given new instance b, its class is l = arg max_{i = 1,⋯, k}x_{ i }. It selects the maximum coefficient in the coefficient vector, and then assigns the class label of the corresponding training instance to this new instance. Essentially, this rule is equivalent to applying nearest neighbor (NN) classifier in the column space of the training instances. In this space, the representations of the training instances are identity matrix. The NN rule can be further generalized to the weighted KNN rule. Suppose a Klength vector $\overline{x}$ accommodates the Klargest coefficients from x, and $\overline{c}$ has the corresponding K class labels. The class label of b can be designated as l = arg max_{i = 1, ⋯, C}s_{i} where ${s}_{i}=\text{sum}\left({\delta}_{i}(\overline{x}\right))$. ${\delta}_{i}\left(\overline{x}\right)$ is a Klength vector and is defined as
The maximum value of K can be k, the number of dictionary atoms. In this case, K is in fact the number of all nonzeros in x. Alternatively, the nearest subspace (NS) rule, proposed in [19], can be used to interpret the sparse coding. NS rule takes the advantage of the discrimination of property in the sparse coefficients. It assigns the class with the minimum regression residual to b. Mathematically, it is expressed as j = min_{1≤i≤C}r_{ i }(b) where r_{ i } (b) is the regression residual corresponding to the ith class and is computed as ${r}_{i}\left(b\right)\mathsf{\text{}}=\phantom{\rule{2.77695pt}{0ex}}{bA{\delta}_{i}\left(x\right)}_{2}^{2}$, where δ_{ i }(x) is defined analogically as in Equation (21).
Algorithm 1 Sparsecodingbased classification
Input: A_{ m×n }: n training instances, c: class labels, B_{ m×p }: p new instances
Output: p: predicted class labels of the p new instances

1.
Normalize each instance to have unit l_{2}norm.

2.
Learn the sparse coefficient matrix X, of the new instances by solving Equation (10), (17), or (19).

3.
Use a sparse interpreter to predict the class labels of new instances, e.g. the NN, KNN, or NS rule.
Optimization
Activeset algorithm for l_{1} LS
The problem in Equation (10) is equivalent to the following nonsmooth unconstrained quadratic programming (QP):
where H_{ k×k }= A^{T}A, and g=  A^{T}b. We thus know that the l_{1}LS problem is a l_{1}QP problem. This can be converted to the following smooth constrained QP problem:
where u is an auxiliary vector variable to squeeze x towards zero. It can be further written into the standard form:
where I is an identity matrix. Obviously, the Hessian in this problem is positive semidefinite as we always suppose H is positive semidefinite in this paper.
A general activeset algorithm for constrained QP is provided in [20], where the main idea is that a working set is updated iteratively until it meets the true active set. In each iteration, a new solution x_{t} to the QP constrained only by the current working set is obtained. If the update step p_{t} = x_{t}  x_{t 1}is zero, then Lagrangian multipliers of the current active inequalities are computed. If all these multipliers corresponding to the working set are nonnegative, then the algorithm terminates with an optimal solution. Otherwise, an active inequality is dropped from the current working set. If the update step p_{t} is nonzero, then an update length α is computed using the inequality of the current passive set. The new solution is updated as x_{ t }= x_{t 1}+ α p_{ t }. If α < 1, then a blocking inequality is added to the working set.
To solve our specific problem efficiently in Equation (24), we have to modify the general method, because i) our constraint is sparse, for the ith constraint, we have x_{ i }  u_{ i } ≤ 0 (if i ≤ k) or x_{ i }  u_{ i } ≤ 0 (if i > k); and ii) when u_{ i }is not constrained in the current working set, the QP constrained by the working set is unbounded, therefore it is not necessary to solve this problem to obtain p_{ t }. In the later situation, p_{t} is unbounded. This could cause some issues in numerical computation. Solving the unbounded problem is timeconsuming if the algorithm is unaware of the unbounded issue. If p_{ t }contains positive or negative ∞, then the algorithm may crash.
We propose the revised activeset algorithm in Algorithm 2 for l_{1}LS sparse coding. To address the potential issues above, we have the following four modifications. First, we require that the working set is complete. That is all the variables in u must be constrained when computing the current update step. (And therefore all variables in x are also constrained due to the specific structure of the constraints in our problem.) For example, if k = 3, a working set {1, 2, 6} is complete as all variables, x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}, are constrained, while {1, 2, 4} is not complete, as u_{3} (and x_{3}) is not constrained. Second, the update step of the variables that are constrained once in the working set are computed by solving the equality constrained QP. The variables constrained twice are directly set to zeros. In the example above, suppose the current working set is {1, 2, 4, 6}, then x_{2}, x_{3}, u_{2}, u_{3} are computed by the constrained QP, while x_{1} and u_{1} are zeros. This is because the only value satisfying the constraint u_{1} = x_{1} = u_{1} is x_{1} = u_{1} = 0. Third, in this example, we do not need to solve the equality constrained QP with four variables. In fact we only need two variables by setting u_{2} = x_{2} and u_{3} = x_{3}. Forth, once a constraint is dropped from the working set and it becomes incomplete, other inequalities must be immediately added to it until it is complete. In the initialization of Algorithm 2, we can alternatively initialize x by 0's. This is much efficient than x= (H)^{1}( g) for largescale sparse coding and very sparse problems.
Activeset algorithm for NNLS and l_{ 1 } NNLS
Both the NNLS problem in Equation (17) and the l_{1}NNLS problem in Equation (19) can be easily reformulated to the following nonnegative QP (NNQP) problem:
Algorithm 2 Activeset l _{1} QP algorithm
Input: Hessian H_{ k×k }, vector g_{k× 1}, scalar λ
Output: vector x which is a solution to min $\frac{1}{2}{x}^{\mathsf{\text{T}}}Hx+{g}^{\mathsf{\text{T}}}x+{\lambda}^{\mathsf{\text{T}}}u,\mathsf{\text{s}}.\mathsf{\text{t}}.u\le x\le u$
% initialize the algorithm by a feasible solution and complete working set
x = (H)^{1}( g); u=  x ;
$\mathcal{R}=\left\{i,\mathsf{\text{}}j\forall i:\mathsf{\text{if}}\phantom{\rule{2.77695pt}{0ex}}{x}_{i}>0\mathsf{\text{let}}\phantom{\rule{0.3em}{0ex}}j=k+i\phantom{\rule{0.3em}{0ex}}\mathsf{\text{otherwise}}j=i\right\}$; % initialize working set
$\mathcal{P}=\left\{\mathsf{\text{1}}:\mathsf{\text{2}}k\right\}\mathcal{R}$; % initialize inactive(passive) set
while true do
% compute update step
let ${\mathcal{R}}_{\mathsf{\text{once}}}$ be the indices of variables constrained once by $\text{}\mathcal{R}$;
p_{2k× 1}= 0;
${p}_{x,{\mathcal{R}}_{\text{once}}}=\mathsf{\text{arg}}\phantom{\rule{0.3em}{0ex}}{\mathsf{\text{min}}}_{q}{q}^{\mathsf{\text{T}}}{H}_{{\mathcal{R}}_{\mathsf{\text{once}}}}q+{\left[{H}_{{\mathcal{R}}_{\mathsf{\text{once}}}}{{x}_{\mathcal{R}}}_{{}_{\mathsf{\text{once}}}}+{g}_{{\mathcal{R}}_{\mathsf{\text{once}}}}+\lambda e\right]}^{\mathsf{\text{T}}}q$ where e_{ i } = 1 if ${u}_{{\mathcal{R}}_{\mathsf{\text{once}}},i}={x}_{{\mathcal{R}}_{\mathsf{\text{once}}},i}$, or 1 if ${u}_{{\mathcal{R}}_{\mathsf{\text{once}}},i}={x}_{{\mathcal{R}}_{\mathsf{\text{once}}},i}$;
if p= 0 then
obtain Lagrange multiplier µ by solving
where A is the constraint matrix in Equation (24)
if ${\mu}_{i}\ge 0\forall i\in \mathcal{R}$then
terminate successfully;
else
$\mathcal{R}=\mathcal{R}j;\mathcal{P}=\mathcal{P}+j\phantom{\rule{2.77695pt}{0ex}}$ where $\phantom{\rule{2.77695pt}{0ex}}j=\mathsf{\text{argmi}}{\mathsf{\text{n}}}_{l\in \mathcal{R}}{\mu}_{l}$;
add other passive constraints to $\text{}\mathcal{R}$ until it is complete;
end if
end if
if p≠ 0 then
$\alpha =\mathsf{\text{min}}\left(1,{\mathsf{\text{min}}}_{i\in \mathcal{P},{a}_{i}^{\text{T}}p\ge 0}\frac{{a}_{i}^{\text{T}}\left[x;u\right]}{{a}_{i}^{\text{T}}p}\right);$
[x; u] = [x; u] + α p;
if α < 1 then
$\mathcal{R}=\mathcal{R}+i$; $\mathcal{P}=\mathcal{P}\phantom{\rule{0.3em}{0ex}}i$. where i corresponds to α;
end if
end if
end while
where H= A^{T}A, g=  A^{T}b for NNLS, and g=  A^{T}b+ λ for l_{1}NNLS.
Now, we present the activeset algorithm for NNQP. This problem is easier to solve than l_{1}QP as the scale of Hessian of NNQP is half that of l_{1}QP and the constraint is much simpler. Our algorithm is obtained through generalizing the famous activeset algorithm for NNLS by [21]. The complete algorithm is given in Algorithm 3. The warmstart point is initialized by the solution to the unconstrained QP. As in Algorithm 2, x can be alternatively initialized by 0's. The algorithm keeps adding and dropping constraints in the working set until the true active set is found.
Algorithm 3 Activeset NNQP algorithm
Input: Hessian H_{ k×k }, vector g_{k× 1}
Output: vector x which is a solution to $\mathsf{\text{min}}\frac{1}{2}{x}^{\mathsf{\text{T}}}Hx+{g}^{\mathsf{\text{T}}}x,\mathsf{\text{s}}.\mathsf{\text{t}}.\phantom{\rule{0.3em}{0ex}}x\ge 0$
x = [(H)^{1}( g)]_{+}; % x= [y]_{+}is defined as x_{ i } = y_{ i } if y_{ i } > 0, otherwise x_{ i } = 0
$\mathcal{R}=\left\{i{x}_{i}=\mathsf{\text{}}0\right\}$; % initialize active set
$\mathcal{P}=\left\{i{x}_{i}>0\right\}$; % initialize inactive(passive) set
µ= Hx+ g; % the lagrange multiplier
while R¹ Æ and min_{i ÎR}(µ_{ i }) < e do
% e is a small positive numerical tolerance
j = arg min_{i ÎR}(µ_{ i }); % get the minimal negative multiplier
$\mathcal{P}=\mathcal{P}+\left\{j\right\};\mathcal{R}=\mathcal{R}\left\{j\right\};$
${t}_{\mathcal{P}}={\left({H}_{\mathcal{P}}\right)}^{1}\left({g}_{\mathcal{P}}\right);$
${t}_{\mathcal{R}}=0;$
while min ${t}_{\mathcal{P}}\le 0$do
$\alpha =\underset{i\in \mathcal{P},{t}_{i}\le 0}{\mathsf{\text{min}}}\frac{{x}_{i}}{{x}_{i}{t}_{i}};$
$\mathcal{K}=\mathsf{\text{arg}}\phantom{\rule{0.3em}{0ex}}{\mathsf{\text{min}}}_{i\in \mathcal{P},{t}_{i}\le 0}\frac{{x}_{i}}{{x}_{i}{t}_{i}}$; % there is one or several indices correspond to α
x = x + α(t  x);
$\mathcal{P}=\mathcal{P}\mathcal{K};\mathcal{R}=\mathcal{R}+\mathcal{K};$
${t}_{\mathcal{P}}={\left({H}_{\mathcal{P}}\right)}^{1}\left({g}_{\mathcal{P}}\right);$
${t}_{\mathcal{R}}=0;$
end while
x = t;
µ = Hx + g;
end while
Parallel activeset algorithms
The formulations of l_{1}QP and NNQP sparse coding for p new instances are, respectively,
and
If we want to classify multiple new instances, the initial idea in [19] and [11] is to optimize the sparse coding one at a time. The interiorpoint algorithm, proposed in [22], is a fast largescale sparse coding algorithm, and the proximal algorithm in [23] is a fast firstorder method whose advantages have been recently highlighted for nonsmooth problems. If we adapt both algorithms to solve our multiple l_{1}QP in Equation (27) and NNQP in Equation (28), it will be difficult to solve the single problems in parallel and share computations. Therefore, the timecomplexity of the multiple problems will be the summation of that of the individual problems. However, the multiple problems can be much more efficiently solved by activeset algorithms. We adapt both Algorithms 2 and 3 to solve multiple l_{1}QP and NNQP in a parallel fashion. The individual activeset algorithms can be solved in parallel by sharing the computation of matrix inverses (systems of linear equations in essence). At each iteration, single problems having the same active set have the same systems of linear equations to solve. These systems of linear equations can be solved once only. For a large value p, that is largescale multiple problems, activeset algorithms have dramatic computational advantage over interiorpoint [22] and proximal [23] methods unless these methods have a scheme of sharing computations. Additionally, activeset methods are more precise than interiorpoint methods. Interiorpoint methods do not allow ${u}_{i}^{2}={x}_{i}^{2}$ and ${u}_{i}^{2}$ must be always greater than ${x}_{i}^{2}$ due to feasibility. But ${u}_{i}^{2}={x}_{i}^{2}$ is naturally possible when the ith constraint is active. u_{ i } = x_{ i } = 0 is reasonable and possible. Activeset algorithms do allow this situation.
Kernel extensions
As the optimizations of l_{1}QP and NNQP only require inner products between the instances instead of the original data, our activeset algorithms can be naturally extended to solve the kernel sparse coding problem by replacing inner products with kernel matrices. The NS decision rule used in Algorithm 1 also requires only inner products. And the weighted KNN rule only needs the sparse coefficient vector and class information. Therefore, the classification approach in Algorithm 1 can be extended to kernel version. For narrative convenience, we also denote the classification approaches using l_{1}LS, NNLS, and l_{1}NNLS sparse coding as l_{1}LS, NNLS, and l_{1}NNLS, respectively. Prefix "K" is used for kernel versions.
Dictionary learning methods
We pursue our dictionarylearningbased approach for biological data, based on the following two motivations. First, since sparsecodingonly approach is a lazy learning, the optimization can be slow for large training set. Therefore, learning a concise dictionary is more efficient for future realtime applications. Second, dictionary learning may capture hidden key factors which correspond to biological pathways, and the classification performance may hence be improved. In the following, we first give the dictionary learning models using Gaussian prior and uniform prior, respectively. Next, we give the classification method based on dictionary learning. We then address the generic optimization framework of dictionary learning. Finally, we show that the kernel versions of our dictionary learning models and the classification approach can be easily obtained.
Dictionary learning models
Now we give our dictionary learning models using Gaussian prior and uniform prior over the dictionary atoms, respectively. Both priors aims to get rid off the arbitrary scale interchange between dictionary and coefficient. Suppose D_{ m×n }is the data of n training instances, and the dictionary A to be learned has k atoms. If the Gaussian prior in Equation (6) is used on the dictionary atom, our dictionary learning models of l_{1}LS, NNLS, and l_{1}NNLS are expressed as follow, respectively:
and
The strength of the Gaussian prior based dictionary learning is that it is flexible to control the scales of dictionary atoms. However, it has two model parameters, which increase the model selection burden in practice. Alternatively, in order to eliminate the parameter α, we design an uniform prior over the dictionary which is expressed as
where p is a constant. That is the feasible region of the dictionary atoms is a hypersphere centered at origin with unit radius, and all the feasible atoms have equal probability. The corresponding dictionary learning models are given in the following equations, respectively:
and
A generic optimization framework for dictionary learning
We devise blockcoordinatedescentbased algorithms for the optimization of the above six models. The main idea is that, in the next step, Y is fixed, and the inner product A^{T}A, rather than A itself, is updated; in the next step, Y is updated while fixing A^{T}A(a sparse coding procedure). The above procedure is repeated until the termination conditions are satisfied.
Now, we show that A can be analytically obtained. For normal prior over dictionary atoms, the optimization of finding A in Equations (29), (30), and (31) is to solve
Taking the derivative with respect to A and setting it to zero, we have
We hence have
where Y^{‡} = Y^{T}(Y Y^{T} + α I)^{1}. The inner product A^{T}A can thus be updated by
We also can compute A^{T}D by
For the uniform prior as in Equation (32), updating unnormalized A while fixing Y in Equations (33), (34), and (35) is to solve the generalized least squares:
Taking derivative with respect to A and setting it to zero, we have
Algorithm 4 The generic dictionary learning framework
Input: K= D^{T}D, dictionary size k, λ
Output: R= A^{T}A, Y
nitialize Y and R= A^{T}A randomly;
r_{ prev } = Inf ; % previous residual
for i = 1 : maxIter do
update Y by solving the activeset based l_{1}LS, NNLS, or l_{1}NNLS sparse coding algorithms;
if Gassian prior over A then
update R= Y^{‡T}D^{T}DY^{‡};
end if
if uniform prior over A then
update R= Y^{†T}D^{T}DY^{†};
normalize R by $R=R./\sqrt{\mathsf{\text{diag}}\left(R\right)\mathsf{\text{diag}}{\left(R\right)}^{\mathsf{\text{T}}}}$;
end if
if i == maxIter or i mod l = 0 then
% check every l iterations
r_{ cur } = f (A, Y ); % current residual of a dictionary learning model
if r_{ prev }  r_{ cur } ≤ e or r_{ cur } ≤ e then
break;
end if
r_{ prev } = r_{ cur };
end if
end for
where Y^{†} = Y^{T}(Y Y^{T})^{1}. The inner products of R= A^{T}A and A^{T}D are computed similarly as for the Gaussian prior. The normalization of R is straightforward. We have $R=R./\sqrt{\mathsf{\text{diag}}\left(R\right)\mathsf{\text{diag}}{\left(R\right)}^{\mathsf{\text{T}}}}$, where ./ and $\text{}\sqrt{\u2022}$ are elementwise operators. Learning the inner product A^{T}A instead of A has the benefits of dimensionfree computation and kernelization.
Fixing A, Y can be obtained via our activeset algorithms. Recall that the sparse coding only requires the inner products A^{T}A and A^{T}D. As shown above, we find that updating Y only needs its previous value and the inner product between training instances.
Due to the above derivation, we have the framework of solving our dictionary learning models as illustrated in Algorithm 4.
Classification approach based on dictionary learning
Now, we present the dictionarylearningbased classification approach in Algorithm 5. The dictionary learning in the training step should be consistent with the sparse coding in the prediction step. As discussed in the previous section, the sparse coding in the prediction step needs the inner products A^{T}A, B^{T}B and A^{T}B which actually is Y^{‡T}D^{T}B or Y^{‡T}D^{T}B.
Algorithm 5 Dictionarylearningbased classification
Input: D_{ m×n }: n training instances, c the class labels, B_{ m×p }: p new instances, k: dictionary size
Output: p: the predicted class labels of the p new instances
{training step:}
1: Normalize each training instance to have unit l_{2} norm.
2: Learn dictionary inner product A^{T}A and sparse coefficient matrix Y of training instances by Algorithm 4.
3: Train a classifier f (θ) using Y(in the feature space spanned by columns of A).
{prediction step:}
1: Normalize each new instance to have unit l_{2} norm.
2: Obtain the sparse coefficient matrix X of the new instances by solving Equation (27), or (28).
3: Predict the class labels of X using the classifier f (θ) learned in the training phase.
Kernel extensions
For Gaussian dictionary prior, the l_{1}LS based kernel dictionary learning and sparse coding are expressed in the following, respectively:
where f (·) is a mapping function. Equations (30), (31), (33), (34), (35) and their sparse coding models can be kernelized analogically. As we have mentioned that the optimizations of the six dictionary learning models, only involves inner products of instances. Thus, we can easily obtain their kernel extensions by replacing the inner products with kernel matrices. Hereafter, if dictionary learning is employed in sparse representation, then prefix "DL" is used before "l_{1}LS", "NNLS", and "l_{1}NNLS". If kernel function other than the linear kernel is used in dictionary learning, then prefix "KDL" is added before them.
Computational experiments
Two highthroughput biological data, including a microarray gene expression data set and a protein mass spectrometry data set, are used to test the performance of our methods. The microarray data set is a collection of gene expression profiles of breast cancer subtypes [24]. This data set includes 158 tumor samples from five subtypes measured on 13582 genes. The mass spectrometry data set is composed of 332 samples from normal class and prostate cancer class [25]. Each sample has 15154 features, that is the masstocharge ratios. Our experiments are separated into two parts. The performance of sparse coding for direct classification is first investigated with respect to accuracy and running time. Then our dimension reduction techniques using dictionary learning are tested.
Sparse coding for direct classification
When dictionary learning was not involved, the dictionary was "lazily" composed by all the training instances available. In our experiment, the activeset optimization methods for l_{1}LS, NNLS, and l_{1}NNLS were tested. The weighted KNN rule and NS rule, mentioned in Algorithm 1, were compared. We set K in the KNN rule to the number of all training instances, which is an extreme case as opposite to the NN rule. Linear and radial basis function (RBF) kernels were employed. We compared our activeset algorithms with the interiorpoint [22] method and proximal [23] method for l_{1}LS sparse coding (abbreviated by l_{1}LSIP and l_{1}LSPX). Benchmark classifiers, including kNN and SVM using RBF kernel, were compared. We employed fourfold crossvalidation to partition a data set into training sets and test sets. All the classifiers ran on the same training and test splits for fair comparison. We performed 20 runs of crossvalidation and recorded the averages and standard deviations. Line or grid search was used to select the parameters of a classifiers.
The average accuracies of all classifiers with the corresponding standard deviations on both data sets are compared in Figure 1, from which we have four observations. First, the weighted KNN rule obtained comparable accuracies with the NS rule. The advantage of the KNN rule over the NS rule is that the former predicts the class labels based on the sparse coefficient vector solely, while the later has to use the training data to compute regression residuals. Therefore the KNN rule is more efficient and should be preferred. Second, on the Prostate data, the sparse coding method l_{1}LS and Kl_{1}LS achieved the best accuracy. This convinces us that sparse coding based classifiers can be very effective for classifying highthroughput biological data. Third, the nonnegative models including NNLS, l_{1}NNLS and their kernel extensions achieved competitive accuracies with the stateoftheart SVM on both data set. Fourth, the l_{1}LS sparse coding using our activeset algorithm had the same accuracies as that using the interiorpoint algorithm and proximal algorithm on Breast data. But on Prostate data, the proximal method yielded a worse accuracy. This implies that our activeset method converges to the global minima as the interiorpoint method, while performance may be deteriorated by the approximate solution obtained by the proximal method in practice.
The mean running time (in second) of crossvalidation are shown in Figure 2. For better comparison, logarithm of base two was taken on the results. First of all, we can clearly see that the interiorpoint method is very slow for the l_{1}LS sparse coding. Second, our activeset method is more efficient than the proximal method on Breast data. This is because i) activeset methods are usually the fastest ones for quadratic and linear programmes of small and median size; and ii) expensive computations, like solving systems of linear equations, can be shared in the activeset method. Third, NNLS and l_{1}NNLS have the same timecomplexity. This is reasonable, because both can be formulated to NNQP problem. These nonnegative models are much simpler and faster than the nonsmooth l_{1}LS model. Hence, if similar performance can be obtained by l_{1}LS and the nonnegative models in an application, we should give preference to NNLS and l_{1}NNLS.
Dictionarylearning for feature extraction
The performance of various dictionary learning models with linear and RBF kernels were investigated on both Breast and Prostate data sets. The Gaussianprior based and uniformprior based dictionary learning models were also compared. Again, our activeset dictionary learning was compared with the interiorpoint [22] method and proximal method [23]. The semiNMF based on multiplicative update rules [26] is also included in the competition. As in the previous experiment, fourfold crossvalidation was used. All methods ran on the same splits of training and test sets. We performed 20 runs of crossvalidation for reliable comparison. After feature extraction by using dictionary learning on the training set, linear SVM classifier was learned on the reduced training set and used to predict the class labels of test instances.
In Figure 3, we show the mean accuracy and standard deviation of 20 results for each method. First, we can see that the models with Gaussian prior on dictionary atoms obtained similar accuracies as the uniform prior. Second, with the comparison to sparse coding methods on Breast data as given Figure 1, we can see that dictionary learning increases the prediction accuracy. Third, from the comparison of Figures 3 and 1, we find that the dictionary learning based methods  DLNNLS and DLl_{1}NNLS, obtained similar accuracies as the sparse coding methods  NNLS and l 1NNLS. This convinces us that dictionary learning is a promising feature extraction technique for highdimensional biological data. On Prostate data, we can also find that the accuracy obtained by DLl_{1}LS is slightly lower than l_{1}LS. This is may be because the dictionary learning is unsupervised. Fourth, using the model parameters, DLl_{1}LS using activeset algorithm obtained higher accuracy than DLl_{1}LSIP and DLl_{1}LSPX on Prostate data. The accuracy of DLl 1LS is also slightly higher than that of DLl_{1}LSIP on Breast data. Furthermore, the nonnegative DLNNLS yielded the same performance as the wellknown semiNMF, while further corroborates the satisfactory performance of our dictionary learning framework. Finally, the kernel dictionary learning models achieved similar performance as their linear counterparts. We believe that the accuracy could be further improved by a suitably selected kernel.
We compare the mean computing time of all the feature extraction methods in Figure 4. First, we can see that DLl_{1}LS using activeset algorithm is much more efficient than DLl_{1}LSIP, DLl_{1}LSPX, and semiNMF using multiplicative update rules. Second, the nonnegative dictionary learning models are more efficient than the l_{1}regularized models. Therefore as in the sparse coding method, priority should be given to the nonnegative models when attempting to use dictionary learning in an application.
Conclusions
In this study, l_{1}regularized and nonnegative sparse representation models are comprehensively studied for the classification of highdimensional biological data. We give a Bayesian treatment to the models. We prove that the sparse coding and dictionary learning models are in fact equivalent to MAP estimations. We use different priors on the sparse coefficient vector and the dictionary atoms, which lead to various sparse representation models. We propose parallel activeset algorithms to optimize the sparse coding models, and propose a generic framework for dictionary learning. We reveal that the sparse representation models only use inner products of instances. Using this dimensionfree property, we can easily extend these models to kernel versions. With the comparison with existing models for highdimensional data, it is shown that our techniques are very efficient. Furthermore, our approaches obtained comparable or higher accuracies. In order to promote the research of sparse representation in bioinformatics, the MATLAB implementation of the sparse representation methods discussed in this paper can be downloaded at [27].
Our Bayesian treatment may inspire the readers to try other prior distributions in order to design new sparse representation models. It also helps to discover the similarity and difference between sparse representation and other dimension reduction techniques. Our kernel versions can also be used to classify tensor data where an observation is not a vector but a matrix or tensor [28]. They can also be applied in the biomedical text mining and interaction/relational data where only the similarities between instances are known.
We will apply our technique to other highthroughput data, such as microarray epigenomic data, gene copy number profiles, and sequence data. We will impose both sparseinducing prior on dictionary atoms and coefficients. Inspired by Bayesian factor analysis, we will investigate the variable selection methods using sparse dictionary. The sparse dictionary analysis would help us to uncover the biological patterns hidden in the highdimensional biological data. Furthermore, combining Bayesian sparse representation and Bayesian regression leads to Bayesian sparse representation regression model, which is very helpful for designing supervised dictionary learning. Finally we should mention that we are working on a decomposition method for sparse coding which is efficient on largescale biological data where there are at least thousands of samples.
References
 1.
Furey T, Cristianini N, Duffy N, Bednarski D, Schummer M, Haussler D: Support vector machine classification and validation of cancer tissue samples using microarray expression data. Bioinformatics. 2000, 16 (10): 906914. 10.1093/bioinformatics/16.10.906.
 2.
Li Y, Ngom A: The regularized linear models and kernels toolbox in MATLAB. Tech. rep., School of Computer Science. 2013, University of Windsor, Windsor, Ontario, [https://sites.google.com/site/rlmktool]
 3.
Wall M, Rechtsteiner A, Rocha L: Singular value decomposition and principal component analysis. A Practical Approach to Microarray Data Analysis. Edited by: Berrar D, Dubitzky W, Granzow M, Norwell, MA. 2003, Kluwer, 91109.
 4.
Carvalho C, Chang J, Lucas J, Nevins J, Wang Q, West M: Highdimensional sparse factor modeling: applications in gene expression genomics. Journal of the American Statistical Association. 2008, 103 (484): 14381456. 10.1198/016214508000000869.
 5.
Elad M: Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. 2010, New York: Springer
 6.
Bengio Y, Courville A, Vincent P: Representation learning: a review and new perspectives. Arxiv. 2012, 1206.5538v2
 7.
Elad M, Aharon M: Image denoising via learned dictionaries and sparse representation. CVPR, IEEE Computer Society. 2006, Washington DC: IEEE, 895900.
 8.
Li Y, Ngom A: The nonnegative matrix factorization toolbox for biological data mining. BMC Source Code for Biology and Medicine. 2013, 8: 1010.1186/17510473810.
 9.
Lee DD, Seung S: Learning the parts of objects by nonnegative matrix factorization. Nature. 1999, 401: 788791. 10.1038/44565.
 10.
Li Y, Ngom A: Fast sparse representation approaches for the classification of highdimensional biological data. Bioinformatics and Biomedicine (BIBM), 2012 IEEE International Conference on: 47 October 2012. 2012, 16. 10.1109/BIBM.2012.6392688.
 11.
Hang X, Wu FX: Sparse representation for classification of tumors using gene expression data. J. Biomedicine and Biotechnology. 2009, 2009: ID 403689
 12.
Rowe D: Multivariate Bayesian Statistics: Models for Source Separation and Signal Unmixing. 2003, Boca Raton, FL: CRC Press
 13.
West M: Bayesian factor regression models in the "large p, small n" paradigm. Bayesian Statistics. 2003, 7: 723732.
 14.
Bruckstein AM, Donoho DL, Elad M: From sparse solutions of systems of equations to sparse modeling of signals and images. SIAM Review. 2009, 51: 3481. 10.1137/060657704.
 15.
Tibshirani R: Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological). 1996, 58: 267288.
 16.
Yin J, Liu X, Jin Z, Yang W: Kernel sparse representation based classification. Neurocomputing. 2012, 77: 120128. 10.1016/j.neucom.2011.08.018.
 17.
Gao S, Tsang IWH, Chia LT: Kernel sparse representation for image classification and face recognition. ECCV. 2010, Springer, 114.
 18.
Olshausen B, Field D: Sparse coding with an overcomplete basis set: a strategy employed by V1?. Vision Research. 1997, 37 (23): 33113325. 10.1016/S00426989(97)001697.
 19.
Wright J, Yang A, Ganesh A, Sastry SS, Ma Y: Robust face recognition via sparse representation. TPAMI. 2009, 31 (2): 210227.
 20.
Nocedal J, Wright SJ: Numerical Optimization. 2006, New York: Springer, 2
 21.
Lawson CL, Hanson RJ: Solving Least Squares Problems. 1995, Piladelphia: SIAM
 22.
Kim SJ, Koh K, Lustig M, Boyd S, Gorinevsky D: An interiorpoint method for largescale l 1regularized least squares. IEEE J. Selected Topics in Signal Processing. 2007, 1 (4): 606617.
 23.
Jenatton R, Mairal J, Obozinski G, Bach F: Proximal methods for hierarchical sparse coding. JMLR. 2011, 12 (2011): 22972334.
 24.
Hu Z: The molecular portraits of breast tumors are conserved across microarray platforms. BMC Genomics. 2006, 7: 9610.1186/14712164796.
 25.
Petricoin EI: Serum proteomic patterns for detection of prostate cancer. J. National Cancer Institute. 2002, 94 (20): 15761578. 10.1093/jnci/94.20.1576.
 26.
Ding C, Li T, Jordan MI: Convex and seminonnegative matrix factorizations. TPAMI. 2010, 32: 4555.
 27.
The sparse representation toolbox in MATLAB. [https://sites.google.com/site/sparsereptool]
 28.
Li Y, Ngom A: Nonnegative matrix and tensor factorization based classification of clinical microarray gene expression data. Bioinformatics and Biomedicine (BIBM), 2010 IEEE International Conference on: 1821 December 2010. 2010, 438443.
Acknowledgements
This article is based on "Fast Sparse Representation Approaches for the Classification of HighDimensional Biological Data", by Yifeng Li and Alioune Ngom which appeared in Bioinformatics and Biomedicine (BIBM), 2012 IEEE International Conference on. © 2012 IEEE [10]. This research has been partially supported by IEEE CIS Walter Karplus Summer Research Grant 2010, Ontario Graduate Scholarship 20112013, and The Natural Sciences and Engineering Research Council of Canada (NSERC) Grants #RGPIN2281172011.
Declarations
The publication costs for this article were funded by Dr. Alioune Ngom with his Natural Sciences and Engineering Research Council of Canada (NSERC) Grants #RGPIN2281172011.
This article has been published as part of BMC Systems Biology Volume 7 Supplement 4, 2013: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2012: Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/7/S4.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
YL proposed the original idea, did literature survey, implemented the methods, conducted the experiments, and wrote the first draft of this paper. AN supervised the whole research, gave constructive suggestions, and finished the final version of this paper.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Published
DOI
Keywords
 Sparse Representation
 Sparse Code
 Training Instance
 Dictionary Learning
 Sparse Coefficient