Fast semi-supervised discriminant analysis for binary classification of large
data-sets
Joris Taverniera,b,∗, Jaak Simmb , Karl Meerbergena , Joerg Kurt Wegnerc , Hugo Ceulemansc , Yves
Moreaub
a Department

of Computer Science, KU Leuven, Celestijnenlaan 200A, B-3001 Heverlee - Belgium
KU Leuven, Kasteelpark Arenberg 10, bus 2446, B-3001 Heverlee - Belgium
c Janssen Pharmaceutica, Turnhoutseweg 30, B-2340 Beerse - Belgium

arXiv:1709.04794v2 [] 1 Mar 2018

b ESAT-STADIUS,

Abstract
High-dimensional data requires scalable algorithms. We propose and analyze three scalable and related
algorithms for semi-supervised discriminant analysis (SDA). These methods are based on Krylov subspace
methods which exploit the data sparsity and the shift-invariance of Krylov subspaces. In addition, the
problem definition was improved by adding centralization to the semi-supervised setting. The proposed
methods are evaluated on a industry-scale data set from a pharmaceutical company to predict compound
activity on target proteins. The results show that SDA achieves good predictive performance and our
methods only require a few seconds, significantly improving computation time on previous state of the art.
Keywords: Semi-supervised learning, Semi-supervised discriminant analysis, large-scale
2010 MSC: 65F15, 65F50, 68T10
1. Introduction
In fields such as data mining, information retrieval or pattern recognition, dimensionality reduction of high-dimensional data is often used for
classification. This assumes that data can be accurately represented in a lower dimensional manifold.
For instance, principal component analysis (PCA)
[1, 2] is a widely used unsupervised feature extraction method. In addition, linear discriminant analysis (LDA) [2, 3, 4] is a supervised method using
dimensionality reduction to classify data.
For certain data sets, only some sample labels
are known. For certain applications, verifying data
point labels is extremely expensive. These applications can benefit from models that incorporate both
labeled and unlabeled data in a semi-supervised setting [5]. Semi-supervised learning is built on the
∗ Corresponding

author
Email addresses: Joris.Tavernier@cs.kuleuven.be
(Joris Tavernier), jaak.simm@esat.kuleuven.be (Jaak
Simm), Karl.Meerbergen@cs.kuleuven.be (Karl
Meerbergen), jwegner@its.jnj.com (Joerg Kurt Wegner),
hceulema@its.jnj.com (Hugo Ceulemans),
moreau@esat.kuleuven.be (Yves Moreau)
Preprint submitted to Elsevier

assumption that taking advantage of the inherent
structure of unlabeled and labeled data aids in a
classification task.
The difficulty and cost of labeling samples for
certain applications stimulated surging interest in
semi-supervised learning. This interest has led to
several different methods of incorporating information of the unlabeled samples to aid the supervised
classification task such as Semi-supervised Orthogonal Discriminant Analysis via label propagation
(SODA) [6], SEmi-supervised Local Fisher discriminant analysis (SELF) [7], Trace Ratio based Flexible Semi-supervised Discriminant Analysis, [8, 9],
Flexible Manifold Embedding (FME)[10] and Semisupervised Linear Discriminant Analysis [11]. Most
of these methods try to incorporate geometric information of the unlabeled data and potentially estimate the labels to incorporate the unlabeled data
in a known supervised method.
Target prediction in chemogenomics is an example for which the semi-supervised setting could
prove beneficial. Chemical compounds are modeled by combinations of molecules. Depending on
the specific target protein, the compound can be
March 2, 2018

active or inactive1 . However, due to existing drugs
and research, only a few activity pairs are known.
Labeling these activity pairs in vitro is expensive.
Predicting which molecule affects a protein aids in
the development of new drugs [12, 13]. This is a
binary classification task. A possible feature set
of a compound are the binary substructure fingerprints, which are the substructures occurring in a
compound [14]. The resulting data matrix is a highdimensional sparse binary matrix.
Scalable algorithms and high performance computing are required for such high-dimensional data.
Algorithms based on QR-decomposition or singular value decomposition have cubic computational
complexity and are unfeasible. This eliminates several state-of-the-art semi-supervised learning methods mentioned before due to computational workload. LDA has been studied for large-scale data
and scalable methods are available in the literature [15, 16]. More specifically, Spectral regression
discriminant analysis (SRDA) [15] solves LDA efficiently and is a special case of spectral regression (SR) [17] consisting of two phases: spectral
analysis (SA) of an embedded graph and regression. Using manifold regularization [18], LDA was
extended to the semi-supervised setting leading to
semi-supervised discriminant analysis (SDA) [19].
Here, we present data centralization to improve
the problem definition of SDA and how to efficiently
implement the centralization. Then, we use the
shift-invariance property of Krylov subspace methods to find solutions for different regularization parameters in a computationally feasible manner. We
suggest three scalable algorithms for SDA and compare predictive performance and execution time.
These algorithms are applied on real-life large-scale
data from the database of Janssen Pharmaceutica.
Our main contributions are

• applying semi-supervised learning in the form
of SDA to target prediction in chemogenomics
and comparing different algorithms for solving
SDA.
Section 2.1 presents a brief review of LDA and
its properties and Section 2.3 builds on the properties of LDA to introduce SDA. In Section 3 the
necessary tools are explained to apply SDA on large
data sets and our scalable algorithms are defined in
Section 4. Finally SDA is applied on two data sets
from chemogenomics in Section 5.
2. Discriminant analysis
In this section, we review LDA and present an existing scalable algorithm for LDA. Next, we detail
a semi-supervised extension for LDA, named SDA.
Finally, centralization is detailed for LDA and provided for SDA.
Symbol
D
N
l
C
Nc
xi
X
SB
SW
ST
L
z
w
SR
SA
CSR
SDA

• an efficient implementation of SDA using the
shift invariance property of Krylov subspace
methods for multiple regularization parameters and additionally improve the SDA problem definition by adding centralization;

Description
the number of features
the total number of data samples
the number of labeled data samples
the number of classes
the number of data samples in class c
i-th data sample
the data matrix X ∈ RN ×D
the between-class scatter matrix
the within-class scatter matrix
the total scatter matrix
the laplacian matrix
eigenvector of the spectral analysis
eigenvalue problem, z ∈ RN ×1
eigenvector of the spectral analysis
eigenvalue problem projected onto the
range of X, w ∈ RD×1
spectral regression
spectral analysis
centralized spectral regression
semi-supervised discriminant analysis

Table 1: Notations and descriptions used in the paper.

2.1. Linear discriminant analysis
LDA [2, 3, 4, 15, 19] is a supervised dimensionality reducing method and finds the directions that
maximally separate the different classes while minimizing the spread within one class. Suppose we
have C classes and N samples xi of dimension D.
The data matrix X ∈ RN ×D is defined as

• applying spectral regression with centralization to SDA and proposing classification of
data using spectral analysis as a standalone
method in the semi-supervised setting;
1 In general, the chemical compound or one of the
molecules the compound exists of can be active on the protein or not.

X = [x1 , . . . , xi , . . . , xN ]T
2

where each row corresponds to a sample. We define
SB as the between-class-scatter matrix
SB =

C
X

2.2. Spectral regression
LDA has been investigated for large-scale data.
Here we present one state-of-the-art algorithm for
LDA for large scale data named Spectral regression
discriminant analysis (SRDA) [15]. SRDA can actually be interpreted as a part of a graph embedding
analysis technique called spectral regression [17].

Nc (µc − µ)T (µc − µ)

c=1

with Nc the number of data points in class c, µ the
mean sample and µc the mean sample of class c.
The within-class-scatter matrix SW is defined by
SW =

Nc
C X
X

2.2.1. Spectral regression discriminant analysis
The LDA eigenvalue problem (2) can be interpreted as the projection of the eigenvalue problem

T

(xn − µc ) (xn − µc ).

c=1 n=1

The LDA objective function is then

W zi = λ i zi

wT SB w
.
max T
w∈RN w SW w

on the range of X [15]. This eigenvalue problem is
quite simple due to the block diagonal structure of
W where values in one block are identical. There
are two eigenvalues: 1 with multiplicity C and 0
with multiplicity N − C[15]. In general, the space
spanned by the eigenvectors with eigenvalue 1 is the
same as the space spanned by the vectors

Using the total scatter matrix ST
ST =

N
X

(xn − µ)T (xn − µ)

n=1

= SB + SW
the objective function can then equivalently be expressed as
wT SB w
.
max T
w∈RN w ST w
Since SB and ST are positive semi-definite, the optimization problem can be reformulated as finding
the large eigenvalues and eigenvectors of
SB wi = λi (ST + βI)wi

zc = ( 0, . . . , 0 , 1, . . . , 1, 0, . . . , 0 )T
| {z } | {z } | {z }
Pc−1

m=1

Nm

Nc

PC

m=c+1

(4)

Nm

due to the special block diagonal structure of W .
In order to gain some insight in eigenvalue problem (3), we will consider the binary case. If we only
have two classes, then the matrix W has rank 2.
Taking a linear combination of z1 and z2 in (4), two
eigenvectors with eigenvalue 1 of eigenvalue problem (3) are

(1)

which we call the LDA eigenvalue problem. The
regularization parameter β is introduced due to
possibility of singularity of the total scatter matrix and as a countermeasure against overfitting.
Suppose the data matrix X is centered and sorted
according to the class of the data points.

znd

X = [X1 ; X2 ; . . . ; XC ]
where Xi are the data samples of class i. The
between-class covariance matrix can then be simplified to
SB = X T W X

 
 
1
1
 .. 
 .. 
.
 . 
 
 
1
1



=   , zd = 
−1
1
 
.
 . 
 .. 
 .. 
1
−1

where zd has value 1 for the samples of class 1 and
value −1 for the samples of class 2. It is obvious
that eigenvector znd with all values equal to 1 has
no discriminative information. The eigenvector wd ,
with Xwd = zd , is the direction that maximizes
the LDA objective function when zd ∈ Range(X).
In general, the eigenvalue associated with wd is
smaller but close to 1 since it is the projection onto
the range of X of the eigenvector zd with eigenvalue
1.

with W an N ×N block diagonal matrix with C matrices (W1 , . . . , WC ) on its diagonal and each value
of matrix Wc ∈ RNc ×Nc is equal to 1/Nc and Nc is
the number of samples in class c. The total scatter
matrix can be written as ST = X T X and the LDA
eigenvalue problem is then
X T W Xwi = λi (X T X + βI)wi .

(3)

(2)
3

SRDA uses the knowledge that the space spanned
by the eigenvectors zc (3) includes the nondiscriminative eigenvector znd . SRDA orthogonalizes the vectors zc with respect to znd , resulting
in C − 1 discriminative eigenvectors z1 , . . . , zC−1
of eigenvalue problem (3). SRDA then solves the
linear systems Xwc = zc for c = 1, . . . , C − 1 in
a regularized least-squares (RLS) manner, resulting in the C − 1 eigenvectors of the LDA eigenvalue
problem which are the directions that maximize the
LDA objective function.

2.3. Semi-supervised discriminant analysis
In the semi-supervised setting, the class of each
data point is not always known and we will refer
to the data points where the class is known as labeled data and unlabeled data otherwise. For both
the unlabeled and labeled data, their structure in
the feature space, xi , is known. Semi-supervised
learning combines both the class information of the
labeled data and the structure in the feature space
of labeled and unlabeled data. SDA [19] uses manifold regularization to extend LDA to the semisupervised setting. Information of the geometric
proximity of the data points in the feature space is
contained in the regularizer.
Using a binary similarity matrix S ∈ RN ×N ,
where sij = 1 = sji if two samples are similar and 0
otherwise, a natural regularizer R(w) for LDA can
be defined as

2.2.2. Spectral regression
SRDA fits in the spectral regression (SR) framework [17]. SR provides a graph embedding interpretation to different dimensionality reduction
techniques such as Semantic Subspace Projection
[20] or Locality Preserving Projections [21]. Based
on a graph, a symmetric matrix W represents the
weights of the edges between two nodes or in our
case data samples. Note that for LDA the resulting
graph consists of the data points and two points are
connected by an edge with weight 1/Nc if they have
the same class c.
Spectral analysis (SA) of the graph computes the
eigenvectors associated with the largest eigenvalues
of
W z = λDz
(5)
PN
with D a diagonal matrix and Dii = n=1 Win . As
shown by Cai et al [17], this corresponds to penalizing two connected nodes i, j with a large difference
|zi − zj |. For new data samples a linear function is
chosen such that zi = xi w, resulting in z = Xw.
Eigenvalue problem (5) then results in
T

T

X W Xw = λX DXw

R(w) =

n
n X
X

sij (wT xi − wT xj )2

i=1 j=1

=2

n
X


(wT xi )2 di −

i=1
T

n
X


sij wT xi wT xj 

j=1
T

= 2w X (D − S)Xw
where di is the sum of the ith row of S and D ∈
Rn×n the diagonal matrix with di on its diagonal.
The matrix L = D −S is the Laplacian matrix of S.
Suppose the number of labeled data is `, the SDA
eigenvalue problem is then


W 0
T
wi =
X
X
0 0
{z
}
|
Between-labeled-class-scatter


I
λi [X (1 − α) `
0
|
{z

(6)

T

which for LDA is similar to (2) without regularization.
SR thus first computes the eigenvectors in (5)
using the Lanczos method [22] and is called the
spectral analysis phase. Note that for SRDA these
eigenvectors were computed analytically and no
eigenvalue solve is required. SR secondly solves the
linear system z = Xw in a regularized least-squares
manner
w = (X T X + βI)−1 X T z


0
X+
0
}

Total-labeled-scatter
T
X
αLX}
| {z

+

βI
|{z}

] wi

(7)

Manifold regularization

Regularization

where the first ` rows of X are the labeled data
points and W ∈ R`×` the same matrix as for LDA.
The regularization parameter α is limited to the interval [0, 1] and resembles how much weight is given
to the unsupervised information. Note that the information of the unlabeled data is only used in the

to find the desired directions w in (2) or (6) during
the regression phase.
4

Laplacian matrix L. SDA can be interpreted under the spectral regression framework and its SA
formulation (5) is






W 0
I` 0
z = λ (1 − α)
+ αL z (8)
0 0
0 0

Note that due to the structure of W ,
X T W Xwnd

= X T W znd
= X T znd
= Nµ

with µ = X T znd /N the mean sample in the feature space. If the data is centralized, the nondiscriminative eigenvector wnd is moved to the common null space of SB and ST of the LDA eigenvalue
problem (1).

and is followed by the usual regression phase.
2.4. Similarity matrix
Semi-supervised learning is based on the assumption that similar data points in the feature space
are more likely to have the same class. Defining
the similarity between data points is an important
part of the semi-supervised setting. SDA uses a
binary similarity matrix but this is not a necessary condition. For a binary similarity matrix, the
first approach we considered is creating a k-nearest
neighbor graph [19, 20].
The value of sij ∈ S is 1 if there exists an undirected edge in the k-nearest neighbor graph from
data point i to data point j. The k-nearest neighbor graph is built using a distance measure in the
feature space. The distance for each data point is
calculated to all the other data points. The data
point is then connected by an undirected edge to
all the other points maximally separated by the distance of the k-nearest neighbor.
The second approach is to create a binary similarity matrix by thresholding. Using a distance measure a threshold is defined for all data points and
if the distance from data point i to data point j is
smaller than the threshold, the data points are similar and sij = sji = 1. Note that using a threshold
instead of the k-nearest neighbors graph may lead
to samples having no neighbors at all.

2.5.2. Centralization for SDA
The centralization trick from LDA to remove the
non-discriminative eigenvector can be applied for
SDA as well. Consider the vector znd ∈ RN ×1
with all values equal to 1. From the definition of
the Laplacian matrix L, we know that Lznd = 0.
Therefore, znd is an eigenvector of the SA eigenvalue problem for SDA






W 0
I 0
znd = λnd (1 − α) `
znd + αLznd
0 0
0 0
(9)
with λnd = 1/(1−α). The SDA eigenvalue problem
can be interpreted as a projection of this eigenvalue
problem onto the range of X. As in LDA, we have


W 0
T
X
z = `µ`
0 0 nd
with µ` the mean of the labeled data in the feature
space. To remove the non-discriminative eigenvector, we centralize the data around the mean of the
labeled data since the Laplacian does not influence
the eigenvector wnd , found by solving the system
Xwnd = znd . This is an improvement to the original formulation of SDA [19].

2.5. Centralization
Centralization, as for PCA, plays an important
role in LDA. We will first look into the details of
the influence of centralization for LDA and extend
the same idea to SDA.

3. Large-scale data
Exact solutions for high dimensional data are
costly to calculate and are no longer feasible. First
we look at an iterative method to solve a large generalized singular eigenvalue problem in Section 3.1
and next we look at Krylov subspace methods to
solve linear systems in Section 3.2.

2.5.1. Centralization for LDA
The non-discriminative eigenvector znd can be
removed from the LDA eigenvalue problem (2)
by centralization of the data [15, 16]. If znd ∈
Range(X), then the corresponding vector wnd , with
Xwnd = znd , is an eigenvector of the LDA eigenvalue problem

3.1. Subspace iteration
For the SDA eigenvalue problem, we have to solve
a singular generalized eigenvalue problem

X T W Xwnd = λi (X T Xwnd + βwnd ).

Av = λBv
5

(10)

with

where A is of low rank (C − 1). It is known that
subspace iteration has geometric asymptotic convergence for eigenvalue λi with ratio λj+1 /λi for
i = 1, . . . , j with j the dimension of the subspace
[23]. Choosing j = C − 1 leads to instant convergence due to the low rank of A (λj+1 = 0). The
subspace iteration algorithm applied to eigenvalue
problem (10) is given in Algorithm 1 with R randomly generated from a uniform distribution between −1 and 1.



γ1
..
.

..

..

..

.

.

.
γk−1





γk−1 
αk

a tridiagonal matrix.
The approximate solution of the linear system
Bw = b is found by searching the space spanned
by the vectors of the Krylov space K. One way to
find an approximate solution is to make the residual
orthogonal to the Krylov base ((B + βI)wk − b ⊥
K). This is true when wk = Kzk where zk is the
solution of the linear system

Algorithm 1 Subspace iteration for (10)
Initialize R ∈ RD×C−1 randomly
Solve V from BV = AR
3: return V
1:
2:

(T + βI)zk = ||b||2 e1

(12)

of dimension k. If the matrix is positive-definite,
the solution minimizes the error under the B-norm.
This method is the Lanczos algorithm to solve a
linear system and is related to the conjugate gradient method (CG, Algorithm 2) for which they use
a Cholesky factorization of T to define a recursive
process [22].

The computational bottleneck is solving the linear system with B. When B is high-dimensional
this can no longer be done exactly in a feasible manner and iterative methods are used. The convergence of subspace iteration is no longer guaranteed
in one iteration. In general, many applications such
as SDA require only an approximate solution and
for these applications it is not necessary to solve the
problem to machine precision.

3.2.2. Regularization
Singular eigenvalue problems are well understood
for small scale matrices [25, 26]. However, these
algorithms are unfeasible for large scale problems,
since singular subspace deflation, removing the null
space from both matrices, is required and has cubic complexity. In contrast, adding a regularization parameter to matrix B (10) is computationally cost-effective, equally well-defined and offers a
countermeasure against overfitting.
Replacing B by B + βI makes the eigenvalue
problem well defined for which we can use a standard method, such as subspace iteration. If β is
small, the eigenvalues and eigenvectors of interest
do not change significantly [27]. In addition, the
singularity of B is inexact, namely, the singular
values quickly decay, but they do not reach zero.
These singular values are the result of ‘noise’ in the
data or a mismatch in the imposed classification.
The effect of this ‘noise’ is reduced by adding the
regularization parameter β.
In general, the best solution is chosen from a
range of β values in the terms of a performance
measure. Krylov bases are invariant to shifts or
regularization parameters and exploiting this property is very efficient [28, 29, 30]. In the previous
discussion, the Krylov vectors are independent of

3.2. Linear systems
We use Krylov subspace based methods for solving a linear system Bw = b on line 2 of Algorithm
1. A Krylov base K is an orthogonal basis determined by the space of a sequence of matrix-vector
multiplications K = span{b, Bb, . . . , B k−1 b}. A
nice property of these methods is that if the matrix
B is sparse the cost of a matrix-vector multiplication is linear with the number of nonzero entries in
the matrix B. For the SDA eigenvalue problem (7),
the matrix B is a product of sparse matrices, which
we do not compute explicitly since this generally
leads to a dense matrix. This subsection is technical and more detailed information can be found in
the literature [22, 24].

3.2.1. Solving the symmetric linear system
For symmetric matrices (which is the case in
SDA) the idea is to find a Krylov base K such that
K∗ BK = T


α1

 γ1
T =



(11)
6

discriminative eigenvector becomes

the value of β and the tridiagonal matrix T (11) is
not influenced by β. This means that the computationally expensive part, that is, building the Krylov
base, can be reused for several β values instead of
restarting the process for each β value. The parameter β only plays a role in solving the small linear
system in (12). The conjugate gradient method for
the shifted linear systems of the form


T 


 !
11T`
I` 0
(IN −
)X
(1 − α)
+ αL X wnd
0 0
`
T 




11T`
I` 0
)X
(1 − α)
+ αL znd
= (IN −
0 0
`

T
11T
= (IN − ` )X (1 − α)z`
`

(B + βI)w = b

1` 1T
)(1 − α)z`
`
`z`
= (1 − α)X T (z` −
)
`
= 0.

= X T (IN −
is given in Algorithm 3 [30, 31]. Superscript β
means that the variable is stored for all different
regularization parameters β in such a way that
scalars (γ) become vectors (γ β ) and vectors (w)
change into matrices (W β ).

Using the centralized multiplication with X T avoids
the non-discriminative eigenvector. This holds for
the centralized multiplication with X as well and
the derivation is similar.

3.2.3. Centralization for a Krylov subspace
In Section 2.5 it was shown that centralization
avoids calculating the non-discriminative eigenvector. Note that if we apply centralization directly
on the matrix X, we will lose the sparsity of the
data matrix and one of the advantages of Krylov
subspace methods. In order to avoid creating a
dense matrix filled with real values, the centralization is applied in the matrix-vector product. Using
znd = 1 ∈ RN ×1 as the vector of all ones and defining z` = 1` ∈ RN ×1 as the vector with value 1
for the labeled samples and 0 otherwise, the mean
labeled sample in the feature space µ` = X T 1` /`
can then be derived. The centralized multiplication
with X can be written as

4. Four algorithms for SDA
Here we present four methods for solving the
SDA eigenvalue problem for binary classification.
First, we propose an efficient implementation of
the spectral regression for SDA (SR-SDA). Next,
we present our improved version of the spectral regression by adding centralization to SR-SDA, resulting in centralized spectral regression for SDA
(CSR-SDA). We then propose spectral analysis as
a standalone method to classify the data (SA-SDA)
in the semi-supervised setting. Finally, we present
our scalable implementation for solving the centralized SDA eigenvalue problem directly (FSDA). The
different algorithms are implemented in C++.
For the evaluation of the algorithms, we use area
under the curve of the receiver operating characteristic (AUC-ROC). In the case of binary classification, there is only one discriminative eigenvector.
To perform the classification using this eigenvector, the data is projected along this vector, namely
Xw = s ∈ RN ×1 , and results in the sample rating
or score. Here we use the rating directly to calculate the AUC-ROC, but other classification methods like support vector machines [32], nearest neighbor [33] or nearest centroid [34, 35] can be used to
classify the data based on the projection.



11T
(IN − ` )X v = (X − 1µT` )v
`
= Xv − 1hµT` , vi
and the centralized multiplication with X T as

T
11T`
(IN −
)X w = (X − 1µT` )T w
`
= X T w − µ` h1T , wi
were the last inner product is the sum of the elements of w.
For the SDA eigenvalue problem, we only have
to implement one centralized multiplication. For
example, using the centralized multiplication with
X T , the matrix-vector product in SDA for the non-

4.1. Spectral regression for SDA
For the spectral regression framework applied
to SDA (SR-SDA), we need to solve the spectral
7

Algorithm 2 Conjugate gradient
1:

2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

Initialize:

Algorithm 3 Shifted conjugate gradient
1:

Initialize:

• r0 = p0 = b

• α0 = 0,W0β = 0, αβ0 = 0

• w=0

• r0 = b, p0 = b, P0β = b

for i = 0, . . . , k do
qi = Bpi
hri ,ri i
γi = − hp
i ,qi i
wi+1 = wi − γi pi
ri+1 = ri + γi qi
if ||ri+1 || < tol then
Break
end if
,ri+1 i
αi+1 = hri+1
hri ,ri i
pi+1 = ri+1 + αi+1 pi
end for

• γ−1 = 1, ζ β−1 = 1, ζ β0 = 1
for i = 0, . . . , k do
qi = Bpi
hri ,ri i
4:
γi = − hp
,q i

2:
3:

i

5:

ζ βi+1

6:

γ βi = γ

=

i

β
ζβ
i−1 ζ i γi−1
β
β
γi αi (ζ i−1 −ζ i )+ζ β
i−1 γi−1 (1−βγi )
ζβ
i+1
i ζβ
i

11:
12:
13:

β
Wi+1
= Wiβ − γ βi Piβ
ri+1 = ri + γi qi
if ||ri+1 || < tol then
Break
end if
,ri+1 i
αi+1 = hri+1
hri ,ri i
pi+1 = ri+1 + αi+1 pi

14:

αβi+1 = αi

7:
8:
9:
10:

15:
16:
17:
18:
19:
20:
21:

β
Pi+1

β
ζβ
i+1 γ i

ζβ
i γi
β
ζ i+1 ri+1

=
+ αβi+1 Piβ
for s = 0, . . . , size(β) do
if ||ζ βi+1 (s)ri+1 || < tol then
Stop updates for β(s)
end if
end for
end for

Figure 1: The Conjugate Gradient algorithm in Algorithm 2 [22] and the Conjugate Gradient algorithm for various shifted
linear systems in Algorithm 3 [30, 31]. The additional operations for shifted CG are highlighted in blue.

8

avoided by looking at the SA phase in SDA (8). The
SA phase is actually an SDA eigenvalue problem (7)
as well with the identity matrix IN as the data matrix X . In order to avoid the non-discriminative
eigenvector znd , we centralize the identity matrix
around the mean of the labeled data (1` /`). Now
we only need to compute one eigenvector and therefore using the power method [23] suffices, which is
basically subspace iteration with only one righthand side. The difference with SR-SDA is that
CSR-SDA computes one eigenvector less during the
spectral analysis phase and therefore requires only
one linear solve during the regression phase. The
full algorithm for centralized spectral regression for
SDA (CSR-SDA) is given is Algorithm 5.

analysis SDA eigenvalue problem (8). The original framework finds all the eigenvectors of the
spectral analysis SDA eigenvalue problem and then
finds the corresponding feature vectors by solving
a linear system in a regularized least squares manner. In order to find the eigenvectors, we use subspace iteration with the block conjugate gradient
method (BCG) [36] as the linear system solver.
Subspace iteration will actually find the eigenspace
Z ∈ RN ×2 containing both the non-discriminative
eigenvector znd and the discriminative eigenvector
zd . SR would normally proceed by computing the
corresponding vectors wnd and wd in the regression
phase. The computed space Z during the spectral
analysis phase consists of linear combinations of znd
and zd . To compare the original SR with our improvements, we need the discriminative eigenvector
to retrieve the rating of the data samples and compute the AUC-ROC.
The discriminative eigenvector can be extracted
by applying Rayleigh-Ritz subspace projection [23].
Using the eigenspace Z, eigenvalue problem (8)
is transformed to a two-by-two eigenvalue problem Z T AZqi = λi Z T BZqi with i = 1, 2 and A,
B shorthands. The discriminative eigenvector is
then the projection of the second largest eigenvector zd = Zq2 . Next we solve the linear system
(X T X + βI)wd = X T zd using shifted CG (Algorithm 3) to find the discriminative eigenvector wd
in the feature space. The full algorithm is given in
Algorithm 4, where we apply the regression phase to
both eigenvectors zd and znd to match the original
algorithm as close as possible. The LSQR method
is advocated for linear systems arising from regularized least squares [37] and used in spectral regression. However, to the best of our knowledge, there
is no specific LSQR variant exploiting the shiftinvariance property of Krylov subspaces 2 . As mentioned, our implementation thus uses shifted CG to
solve the system for multiple regularization parameters.

4.3. Spectral analysis for SDA
In the semi-supervised setting, the data samples
requiring classification are often available as unlabeled data. The regression step in CSR-SDA solves
the linear system z = Xw in a regularized leastsquares manner (Algorithm 5, line 5). The rating
of the samples is next calculated as the projection
of the data point along the discriminative direction
s = Xw. For the rating of the known unlabeled
data the eigenvector z can be used (s = z), making
the regression step obsolete as shown in Algorithm
6. We added regularization to the SA phase as
countermeasure against overfitting. Spectral analysis for SDA (SA-SDA) as standalone method only
works when the samples to be classified are known
and present in X. If α = 0, the Laplacian is neutralized and the unlabeled data receive a rating of
zero. Therefore, this approach only works if α 6= 0.

4.4. Fast semi-supervised discriminant analysis
Since we are working with sparse data and Krylov
subspace methods, the complexity of the algorithms
are proportional to the number of nonzeros. Hence,
we propose an efficient implementation to solve
the centralized SDA eigenvalue problem directly by
using the power method. The full algorithm for
fast semi-supervised discriminant analysis (FSDA)
is given in Algorithm 7.

4.2. Centralized spectral regression for SDA
SR as detailed above computes the redundant
non-discriminative eigenvector znd which can be
2 The spectral regression code [17] of the authors is available at http://www.cad.zju.edu.cn/home/dengcai/Data/
SR.html and is for MATLAB which uses the function eigs.
Solving the spectral analysis eigenvalue problem for SDA
using a toy data set with eigs already took twice as long as
solving the spectral SDA eigenvalue problem and the linear
system with our own C++ code.

4.5. Comparison of the SDA algorithms
All four algorithms solve the same underlying
problem. In a perfect world without regularization
and if the solution of zd = Xwd exists, these algorithms would have the same performance. In a
9

Algorithm 4 SR-SDA
Given α, β
Initialize R ∈ RN ×2 randomly generated from a uniform distribution between −1 and 1
3: Solve Z ∈ RN ×2 with BCG from






I` 0
W 0
(1 − α)
+ αL Z =
R
0 0
0 0
1:
2:

4:

5:
6:

7:
8:

Solve two-by-two eigenvalue problem




W 0
I
ZT
Zqi = λi Z T (1 − α) `
0 0
0
for i = 1, 2 do
Solve wi with shifted CG from



0
+ αL Zqi
0

(X T X + βI)wi = X T (Zqi )

end for
return Rating: s = Xw2

Algorithm 5 CSR-SDA
Given α, β
Initialize r ∈ RN ×1 randomly generated from a uniform distribution between −1 and 1
<1T ,r>
3: Orthogonalize: r = r − 1 ``
4: Solve z with CG from



 



11T
I 0
W 0
+ αL IN z =
(IN − ` )T (1 − α) `
r
0 0
0 0
`
1:
2:

5:

Solve w with shifted CG from
(X T X + βI)w = X T z

6:

return Rating: s = Xw

Algorithm 6 SA-SDA
Given α, β
Initialize r ∈ RN ×1 randomly generated from a
<1T ,r>
3: Orthogonalize: r = r − 1 ``
4: Solve z with shifted CG from



11T
I
(IN − ` )T (1 − α) `
0
`
1:
2:

5:

uniform distribution between −1 and 1





0
W
+ αL IN + βIN z =
0
0

return Rating: s = z

10


0
r
0

Algorithm 7 FSDA
1:
2:
3:
4:

5:

Given α, β
Initialize r ∈ RD×1 randomly generated from a uniform distribution between −1 and
Calculate mean labeled sample: `µ` = X T 1`
Solve w with shifted CG from







I` 0
W
T T
T T
(X − 1µ` ) (1 − α)
+ αL X + βID w = (X − 1µ` )
0 0
0

1


0
Xr
0

return Rating: s = Xw

real world regularization is required and iterative
methods are used to solve linear systems within a
certain tolerance. The difference between SR-SDA
and CSR-SDA is only computational and these two
algorithms result in the same classification. The
differences between CSR-SDA, SA-SDA and FSDA
are more subtle. FSDA solves the SDA eigenvalue
problem (7) directly and the tolerance used in CG
reflects the true error of the eigenvector. For CSRSDA using the same tolerance in both linear systems does not guarantee this tolerance for the SDA
eigenvalue problem.
All presented iterative methods build a Krylov
subspace and each iteration is linear in the number of samples, features and nonzeros of the matrices. Furthermore, centralization costs 4 linear
operations for each iteration with linear the column
size of the matrix used. Table 2 shows the complexity of the different algorithms given in this section.
Applying centralization for CSR-SDA results in a
speed-up of approximately C/(C − 1) with respect
to SR-SDA, reducing execution time by almost a
factor 2 for binary classification.
The table also shows the main difference between
CG and Shifted CG. We did not take into account
the influence of the β-value on the number of iterations required to reach a certain tolerance and
assumed that each solve takes the same number of
iterations. CG has to rebuild the Krylov subspace
for each value of β, while shifted CG reuses the base
resulting in additional linear operations. Note that
the number of nonzeros nz is generally significantly
larger than the column size of the matrix, and thus
shifted CG results in a speed-up approximately the
number of regularization parameters used.

Algorithm

Complexity

SR-SDA
CSR-SDA

k1 C(nzL +10N )+k2 C(2nzX +10D)
k1 (C − 1)(nzL + 14N ) + k2 (C −
1)(2nzX + 10D)
k1 (C − 1)(nzL + 14N )
k2 (C − 1)(2nzX + nzL + 14D)
k1 Nβ (nz + 10N )
k1 (nz + Nβ 10N + Pβ )

SA-SDA
FSDA
CG
Shifted CG
with
k1 : the number of CG-iterations with dimension N
k2 : the number of CG-iterations with dimension D
nzL : the number of nonzeros of L
nzX : the number of nonzeros of X
Nβ : the number regularization parameters β
Pβ : extra cost related to Nβ used in shifted CG
Table 2: The complexity of the 4 SDA algorithms.

the data sets. Since these methods solve the same
underlying problem, we are interested in convergence speed and difference in the algorithms predictive performances.
This work is performed with drug-protein activity prediction in mind. We have tested our methods
on different targets for two different data sets. The
first data set is described in section 5.1 and is publicly available. The implementation for this data set
was compiled with gcc 5.3.1 and OPENMP 4.0 with
compile option -O3 and executed on a machine with
an Intel(R) Core(TM) i7-6560U (2.20GHz) processor with an L3 cache memory of 4096 Kb and 16
GB of DRAM where all 4 cores were used 3 . The
second industry-scale data set is described in section 5.2 and stronger hardware was required. These

5. Experimental results
3 Implementation and the public ChEMBL data set can
be found here: https://people.cs.kuleuven.be/~joris.
tavernier/.

For evaluating the different methods, we present
the AUC-ROC in function of the execution time for
11

89 ± 4%. In addition, to evaluate the convergence
and execution time for the different methods, we
have randomly selected 6 targets with 2000 or more
labeled compounds. The number of active and inactive compounds for each target is given in Table
4.
For all methods, we varied the number of CG iterations which increases the complexity, but also
the accuracy of CG. Note that the spectral regression methods require two linear solves and the best
pair (execution time, AUC-ROC) is reported. The
different number of iterations for the linear solvers
are given in Table 5. 5-Fold nested cross-validation
was used for the experiments, with the inner-folds
used to pick the best value of the regularization parameter β. Starting from β = 10−9 , we let β vary
by multiplying β with 10 until β = 103 . We report the average execution time for the outer crossvalidation using only one optimal β. The reported
AUC-ROC is the average along with the standard
deviation for all outer folds and 5 different seed initializations of the random number generator.
Figure 3 shows the AUC-ROC in function of the
execution time for the different methods. Solving
SDA directly in the feature space (FSDA) gives
higher or equal average AUC-ROC for all 6 targets
in comparison with the spectral regression methods. There is no difference in AUC-ROC between
SR-SDA and CSR-SDA, but Figure 3 clearly shows
that CSR-SDA requires less execution time to reach
the same predictive performance. For this data,
the AUC-ROC in function of execution times is
higher for FSDA than for CSR-SDA. Ideally, we
prefer to reach high AUC-ROC with the smallest
execution time (top left of the graph). We conclude
that FSDA is more suitable than CSR-SDA for the
ChEMBL data set. Note that SA-SDA reaches high
AUC-ROC values in even less execution time since
there is no regression step. Currently, we do not
know if the difference in AUC-ROC after convergence between CSR-SDA, SA-SDA and FSDA is
specifically for the data or inherently present due
to the way the direction w is calculated.
Table 6 shows the execution time for CG and
shifted CG solving the regularized least squares system in CSR-SDA (Algorithm 5, line 5) with a tolerance of 10−3 and several β-values. Using shifted
CG leads to a speed-up of 7.16 in computation time
for 12 β-values. Note that for the larger β-values
the number of iterations of CG required to reach
the desired tolerance is less than for the small βvalues. Higher regularization parameters result in

experiments were run using a computing node with
dual Intel(R) Xeon(R) E5-2699 v3 (2.30GHz) processors with a L3 cache memory of 46080 Kb and
264 gb of DRAM. For these experiments we used
12 cores.
5.1. ChEMBL data
For drug-protein activity prediction, one of the
possible indicators is the IC50 value. IC50 or the
half maximal inhibitory concentration measures the
substance concentration required to inhibit the activity of a protein by 50%. Using the bioactivity
database ChEMBL version 19 [38], we selected protein targets that had at least 200 or more IC50 values for the compounds. A compound was considered active with an IC50 value less than 1000 nM.
We then kept the protein targets that had at least
30 actives and inactives, resulting in 330 targets
with enough labels to perform 5-fold nested crossvalidation.
For the compounds the extended-connectivity
fingerprints (ECFP [14]) were computed using rdkit [39] with 3 layers. In summary, ECFP generates
the features for a compound in layers. In the first
layer each individual atom receives an initial identifier based on its atomic symbol. In the following
layer this value is concatenated with the bound and
the identifiers of the neighboring atoms and then
hashed to a new identifier for this atom. This process is repeated for a fixed number of layers. These
generated substructures represented by their hash
value are the features for this compound. Duplicate identifiers or duplicate substructures are removed from the feature list. Figure 2 depicts the
creation of the ECFP features for a small chemical compound. The characteristics of the resulting
ECFP data matrix X are given in Table 3.
Parameter
Compounds
Features
Nonzeros
Density

Value
167 668
291 714
12 246 376
0.025038%

Table 3: The characteristics of the ChEMBL data matrix.

The similarity matrix was made using a 5-nearest
neighbor graph and the Tanimoto similarity. The
cross-validated mean and standard deviation for
AUC-ROC of FSDA over the 330 protein targets is
12

H

H

H

C

C

O

O

H

H

Layer 1

O

C

O

C

[O,1,C]

[C,1,O,1,C]

Layer 2

O

C

[O,1,C]

[C,1,O,1,C]

[[O,1,C],1,[C,1,C,1,C,1,O]] Layer 3

H

(b) Resulting features

(a) Compound structure

Figure 2: Example of the ECFP feature generation.

Name
CHEMBL100503
CHEMBL10062
CHEMBL100670
CHEMBL1007
CHEMBL100763
CHEMBL100875

Active

Inactive

Known labels

% labeled

2225
1365
1621
1409
817
2514

1729
1275
491
1170
3144
516

3954
2640
2112
2579
3961
3030

2.36
1.57
1.26
1.54
2.36
1.81

Table 4: The characteristics of the chosen targets for the 167668 chemical compounds of the ChEMBL data set.

reach the same predictive performance with a small
advantage for FSDA as seen for the ChEMBL data
set as well.
The relatively large assays for the Janssen data
allows us to subsample the labeled data. Figure 5
gives the AUC-ROC in function of randomly subsampled labeled data for different values of the parameter α. Using only 0.1 to 1% labels, SDA outperforms LDA for all the targets. Incorporating
information of the unlabeled data trough SDA results in larger predictive performance if the number
of labels is small.

a smaller condition number and faster convergence
for CG [22]. Table 6 additionally shows the speedup if β would not vary significantly, resulting in a
larger speed-up of 8.98. If the β-values do not vary
drastically, the number of iterations of CG do not
change significantly and the achieved speed-up is
larger.
5.2. Janssen data
Next, we investigate the performance of the proposed methods on a industrial-scale data set containing a few million of compounds. This largescale data set was provided by Janssen Pharmaceutica. Out of more than 1000 protein targets,
we chose 6 with a large number of IC50 measurements, see Table 7 for an overview of these targets.
For these targets, a compound was considered active with an IC50 value below 1.0 µM . Using ECFP
over 10 million compound features were generated
and the SDA similarity matrix was created using
the Tanimoto similarity with a threshold of 0.4.
Figure 4 shows the AUC-ROC in function of the
execution time for the targets. The number of iterations for the linear solves was set to 70. As can be
seen from the figure, the centralization improves the
execution time by a factor 2 in comparison with SRSDA. SA-SDA even gains a factor 7-8. All methods

6. Discussion
Semi-supervised learning has mainly been developed for cases where there is a small number of
labels. In our application, there is a medium number of labels and an abundance of unlabeled data
samples. Further research is required specifically
for this kind of large-scale data sets in the semisupervised setting. The information of the unlabeled data is incorporated in the classification task.
To include this information, SDA uses geometric
distance in the feature space to define similarity
between samples. Defining similarity in the feature
space is sensitive to noise and redundant information of the defined features.
13

0.94

0.92

0.92

0.9

0.9
AUC-ROC

AUC-ROC

0.94

0.88
0.86

0.82
10−1

0.86
0.84

FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.84

0.88

FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.82
0.8
10−1

100
Time(s)

100
Time(s)

(a) Target CHEMBL100503

(b) Target CHEMBL10062
0.9

0.9

AUC-ROC

AUC-ROC

0.85
0.85

0.8
FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.75

10−1

0.8

FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.75

0.7

100

10−1

100
Time(s)

Time(s)

(c) Target CHEMBL100670

(d) Target CHEMBL1007

0.88
0.86

0.9
0.85

0.82

AUC-ROC

AUC-ROC

0.84

0.8
0.78
FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.76
0.74

0.8
0.75
FSDA (ours)
CSR-SDA (ours)
SA-SDA (ours)
SR-SDA

0.7
0.65

10−1

100

10−1

Time(s)

100
Time(s)

(e) Target CHEMBL100763

(f) Target CHEMBL100875

Figure 3: Comparison of spectral regression methods and FSDA for the 6 chosen targets of the ChEMBL data set. The
AUC-ROC is given in function of the execution time in seconds. Converging to large AUC-ROC in small execution time (top
left of the graph) is preferable.

14

Method

Description

Iterations

FSDA

Iterations linear solve for SDA

[2, 3, 5, 10, 20, 40, 60, 80]

CSR-SDA

Iterations linear solve for spectral SDA
Iterations RLS solve with X (regression)

[2, 3, 5, 10, 20, 40, 60, 80]
[2, 3, 5, 10, 20, 40, 60, 80]

SA-SDA

Iterations linear solve for spectral SDA

[3, 5, 10, 20, 40, 60, 80]

SR-SDA

Iterations linear solve for spectral SDA
Iterations RLS solve with X (regression)

[2, 3, 5, 10, 20, 40, 60, 80]
[2, 3, 5, 10, 20, 40, 60, 80]

Table 5: The different parameters used in the experiments for each method.

Nβ

β

CG

shifted CG

Speed-up

12

[1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4,
1e-3, 1e-2, 1e-1, 1, 10, 1e2 ]
[1.0e-6, 1.1e-6, 1.2e-6, 1.3e-6, 1.4e-6, 1.5e-6,
1.6e-6, 1.7e-6, 1.8e-6, 1.9e-6, 2.0e-6, 2.1e-6]

43.11 s

6.02 s

7.16

50.98 s

5.68 s

8.98

12

Table 6: The computation time of CG versus shifted CG for the regression phase using a tolerance of 10−3 , the ChEMBL data
matrix and given β-values.

Name

Active

Inactive

Total labels

% labeled

389
685
739
1448
1736
1833

20000
10000
5000
2500
5000
2500

275000
475000
750000
300000
500000
300000

295000
485000
755000
302500
505000
302500

15
25
40
15
25
15

advantage to solve SDA efficiently for different regularization parameters β.
We have applied the spectral regression framework to SDA and showed how to centralize the data
in spectral analysis. We additionally showed how to
use SA as a standalone classification method in the
semi-supervised setting. We applied SDA for target prediction in chemogenomics using data with
medium number of labeled data and an abundance
of unlabeled data.
Using our proposed methods and implementations, it was possible to apply SDA on real-life
large-scale data from the pharmaceutical industry
in seconds, improving on previous state of the art
by a factor 2 or more. Solving the problem directly
in the feature space had experimentally for almost
all targets better or equal predictive performance
than the spectral regression methods.

Table 7: The characteristics of the Janssen targets.

For the chemogenomics application, the features
are exact but do contain redundant information.
The sample similarity could be further improved by
removing non-informative features as suggested in
[11]. Note that for our application, even the smallest speed-up for an algorithm is multiplied by the
number of specific targets, resulting in a significant
reduction of the total execution time.

Conflict of interest
7. Conclusion
No conflict of interest.
We have presented semi-supervised discriminant
analysis for large-scale data. We have shown how
to use the centralization to avoid the calculation of
the non-discriminative eigenvector and an efficient
way to implement the centralization. We used the
shift-invariance property of Krylov subspaces to our

Acknowledgement
The authors would like to thank the reviewers for their comments and suggestions. The authors would additionally like to thank Jean-Marc
15

FSDA (ours), 70 iterations

CSR-SDA (ours), 70 iterations

SR-SDA, 70 iterations

SA-SDA (ours), 70 iterations

0.944
0.944

0.940

0.942

0.982

0.940

0.938

AUC-ROC

AUC-ROC

AUC-ROC

0.984

0.942

0.936
0.934

0.978

0.932

0.938

0.980

0.930
0.936

0.976
10

20

30

40

50
60
Time(s)

70

80

90

10

20

30

40

(a) Target 389

50
60
Time(s)

70

80

90

100

10

30

40
50
Time(s)

60

70

80

(c) Target 739

0.915

0.910
0.988

0.910

0.890

0.905

0.986

AUC-ROC

AUC-ROC

0.900
AUC-ROC

20

(b) Target 685

0.984

0.895
0.890

0.982

0.880

0.900

0.885
0.980

0.870
10

20

30

40

50
60
Time(s)

70

80

90

100

10

20

30

(d) Target 1448

40

50
60
Time(s)

(e) Target 1736

70

80

90

100

0.880

10

20

30

40
50
Time(s)

60

70

80

(f) Target 1833

Figure 4: Comparison of spectral regression methods and FSDA for the 6 Janssen targets. The AUC-ROC is given in function
of the execution time in seconds with a fixed number of iterations equal to 70. Our methods improve the execution time by a
factor 2 in comparison with SR-SDA or more for SA-SDA.

Neefs for annotating assays with protein and gene
identifiers and Vladimir Chupakhin for helping
with the compound normalization of the Janssen
data. This work was supported by Flanders Innovation & Entrepreneurship (IWT, No. IWT130406)
and Research Foundation - Flanders (FWO, No.
G079016N).

[9]

[10]

[11]

References
[12]

References
[1] H. Hotelling, Analysis of a complex of statistical variables into principal components., Journal of Educational Psychology 24 (6) (1933) 417–441.
[2] C. M. Bishop, Pattern Recognition and Machine Learning, Springer-Verlag New York, 2006.
[3] G. J. McLachlan, Discriminant Analysis and Statistical
Pattern Recognition, John Wiley & Sons, Inc., 2005.
[4] T. Hastie, R. Tibshirani, J. Friedman, The Elements of
Statistical Learning, Springer-Verlag New York, 2009.
[5] X. Zhu, A. B. Goldberg, Introduction to semisupervised learning, Synthesis lectures on artificial intelligence and machine learning 3 (1) (2009) 1–130.
[6] F. Nie, S. Xiang, Y. Jia, C. Zhang, Semi-supervised
orthogonal discriminant analysis via label propagation,
Pattern Recognition 42 (11) (2009) 2615 – 2627.
[7] M. Sugiyama, T. Idé, S. Nakajima, J. Sese, Semisupervised local fisher discriminant analysis for dimensionality reduction, Machine Learning 78 (1) (2009) 35.
[8] H. Wang, S. Yan, D. Xu, X. Tang, T. Huang, Trace
ratio vs. ratio trace for dimensionality reduction, in:

[13]

[14]

[15]

[16]

[17]

16

2007 IEEE Conference on Computer Vision and Pattern
Recognition, 2007, pp. 1–8.
Y. Huang, D. Xu, F. Nie, Semi-supervised dimension
reduction using trace ratio criterion, IEEE Transactions
on neural networks and learning systems 23 (3) (2012)
519–526.
F. Nie, D. Xu, I. W. H. Tsang, C. Zhang, Flexible manifold embedding: A framework for semi-supervised and
unsupervised dimension reduction, IEEE Transactions
on Image Processing 19 (7) (2010) 1921–1932.
S. Wang, J. Lu, X. Gu, H. Du, J. Yang, Semi-supervised
linear discriminant analysis for dimension reduction and
classification, Pattern Recognition 57 (2016) 179 – 189.
M. J. Keiser, B. L. Roth, B. N. Armbruster, P. Ernsberger, J. J. Irwin, B. K. Shoichet, Relating protein
pharmacology by ligand chemistry, Nature biotechnology 25 (2) (2007) 197–206.
M. J. Keiser, V. Setola, J. J. Irwin, C. Laggner, A. I.
Abbas, S. J. Hufeisen, N. H. Jensen, M. B. Kuijer, R. C.
Matos, T. B. Tran, et al., Predicting new molecular
targets for known drugs, Nature 462 (7270) (2009) 175–
181.
D. Rogers, M. Hahn, Extended-connectivity fingerprints, Journal of Chemical Information and Modeling
50 (5) (2010) 742–754, pMID: 20426451.
D. Cai, X. He, J. Han, Srda: An efficient algorithm for
large-scale discriminant analysis, IEEE Transactions on
Knowledge and Data Engineering 20 (1) (2008) 1–12.
L. Zhu, D.-S. Huang, A rayleigh–ritz style method for
large-scale discriminant analysis, Pattern Recognition
47 (4) (2014) 1698 – 1708.
D. Cai, X. He, J. Han, Spectral regression: A unified
subspace learning framework for content-based image
retrieval, in: Proceedings of the 15th ACM International Conference on Multimedia, MM ’07, ACM, New

0.950

0.950
0.900

0.850
LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.800

10−1

100
% Labeled

AUC-ROC

AUC-ROC

0.900
0.850
0.800

LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.750
0.700
10−1

101

(a) Target 389

101

(b) Target 685

1.000

0.900

0.800

LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.700

10−1

100
% Labeled

AUC-ROC

0.900
AUC-ROC

100
% Labeled

0.800

LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.700

0.600

10−1

101

(c) Target 739

100
% Labeled

101

(d) Target 1448

1.000
0.900
0.950

0.850
LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.800
0.750
0.700
10−1

100
% Labeled

AUC-ROC

AUC-ROC

0.900

0.800

LDA
SDAα=0.1
SDAα=0.2
SDAα=0.3
SDAα=0.4
SDAα=0.5

0.700

0.600

10−1

101

(e) Target 1736

100
% Labeled

101

(f) Target 1833

Figure 5: The AUC-ROC for the Janssen data in function of number of used labels (randomly subsampled) for 6 different
targets. Using only 0.1 to 1% labels, SDA outperforms LDA for all the targets.

17

York, NY, USA, 2007, pp. 403–412.
[18] M. Belkin, P. Niyogi, V. Sindhwani, Manifold regularization: A geometric framework for learning from
labeled and unlabeled examples, Journal of machine
learning research 7 (Nov) (2006) 2399–2434.
[19] D. Cai, X. He, J. Han, Semi-supervised discriminant
analysis, in: 2007 IEEE 11th International Conference
on Computer Vision, 2007, pp. 1–7.
[20] J. Yu, Q. Tian, Learning image manifolds by semantic subspace projection, in: Proceedings of the 14th
ACM International Conference on Multimedia, MM ’06,
ACM, New York, NY, USA, 2006, pp. 297–306.
[21] X. He, P. Niyogi, Locality preserving projections, in:
S. Thrun, L. K. Saul, B. Schölkopf (Eds.), Advances in
Neural Information Processing Systems 16, MIT Press,
2004, pp. 153–160.
URL
http://papers.nips.cc/paper/
2359-locality-preserving-projections.pdf
[22] J. Demmel, Applied Numerical Linear Algebra, Society
for Industrial and Applied Mathematics, 1997.
[23] B. Parlett, The Symmetric Eigenvalue Problem, Society
for Industrial and Applied Mathematics, 1998.
[24] L. N. Trefethen, D. Bau III, Numerical Linear Algebra,
Vol. 50, Siam, 1997.
[25] J. Demmel, B. Kågström, The generalized schur decomposition of an arbitrary pencil A − λB—robust software
with error bounds and applications. part i: Theory and
algorithms, ACM Trans. Math. Softw. 19 (2) (1993)
160–174.
[26] J. Demmel, B. Kågström, The generalized schur decomposition of an arbitrary pencil A − λB — robust software with error bounds and applications. part ii: Software and applications, ACM Trans. Math. Softw. 19 (2)
(1993) 175–201.
[27] A. Tikhonov, V. Y. Arsenin, Methods for solving illposed problems, John Wiley and Sons, Inc, 1977.
[28] J. van den Eshof, G. L. Sleijpen, Accurate conjugate
gradient methods for families of shifted systems, Applied Numerical Mathematics 49 (1) (2004) 17 – 37,
numerical Algorithms, Parallelism and Applications.
[29] A. Frommer, P. Maass, Fast cg-based methods for
tikhonov–phillips regularization, SIAM Journal on Scientific Computing 20 (5) (1999) 1831–1850.
[30] B. Jegerlehner, Krylov space solvers for shifted linear
systems, arXiv preprint hep-lat/9612014.
[31] A. Sciarra, CG algorithm for shifted matrices.
URL
http://code.compeng.uni-frankfurt.de/
attachments/178/CG-M.pdf
[32] V. N. Vapnik, An overview of statistical learning theory,
IEEE Transactions on Neural Networks 10 (5) (1999)
988–999.
[33] T. Cover, P. Hart, Nearest neighbor pattern classification, IEEE Transactions on Information Theory 13 (1)
(1967) 21–27.
[34] R. M. McIntyre, R. K. Blashfield, A nearest-centroid
technique for evaluating the minimum-variance clustering procedure, Multivariate Behavioral Research 15 (2)
(1980) 225–238.
[35] R. Tibshirani, T. Hastie, B. Narasimhan, G. Chu, Diagnosis of multiple cancer types by shrunken centroids of
gene expression, Proceedings of the National Academy
of Sciences 99 (10) (2002) 6567–6572.
[36] D. P. O’Leary, The block conjugate gradient algorithm
and related methods, Linear Algebra and its Applications 29 (1980) 293 – 322, special Volume Dedicated to

Alson S. Householder.
[37] C. C. Paige, M. A. Saunders, Lsqr: An algorithm for
sparse linear equations and sparse least squares, ACM
Trans. Math. Softw. 8 (1) (1982) 43–71.
[38] A. P. Bento, A. Gaulton, A. Hersey, L. J. Bellis,
J. Chambers, M. Davies, F. A. Krger, Y. Light, L. Mak,
S. McGlinchey, M. Nowotka, G. Papadatos, R. Santos,
J. P. Overington, The chembl bioactivity database: an
update, Nucleic Acids Research 42 (D1) (2014) D1083–
D1090.
[39] RDKit: Open-source cheminformatics, http://www.
rdkit.org, [Online; accessed 18-Januari-2017].

18

