IEICE Transactions on Information and Systems,
vol.E93-D, no.10, pp.2690–2701, 2010.
Revised on June 26, 2011.

1

Superfast-Trainable Multi-Class Probabilistic Classifierby Least-Squares Posterior Fitting
Masashi Sugiyama (sugi@cs.titech.ac.jp)
Tokyo Institute of Technology
andJapan Science and Technology Agency

Abstract
Kernel logistic regression (KLR) is a powerful and ﬂexible classiﬁcation algorithm,
which possesses an ability to provide the conﬁdence of class prediction. However, its training—typically carried out by (quasi-)Newton methods—is rather timeconsuming. In this paper, we propose an alternative probabilistic classiﬁcation
algorithm called Least-Squares Probabilistic Classifier (LSPC). KLR models theclass-posterior probability by the log-linear combination of kernel functions and its
parameters are learned by (regularized) maximum likelihood. In contrast, LSPC
employs the linear combination of kernel functions and its parameters are learned
by regularized least-squares ﬁtting of the true class-posterior probability. Thanks
to this linear regularized least-squares formulation, the solution of LSPC can becomputed analytically just by solving a regularized system of linear equations in aclass-wise manner. Thus LSPC is computationally very eﬃcient and numericallystable. Through experiments, we show that the computation time of LSPC is fasterthan that of KLR by orders of magnitude, with comparable classiﬁcation accuracy.
Keywords
Probabilistic classiﬁcation, kernel logistic regression, class-posterior probability, squaredloss.

1

Introduction

The support vector machine (SVM) [7, 33] is a popular method for classiﬁcation. Various
computationally eﬃcient algorithms for training SVM with massive datasets have been
proposed so far (see [24, 16, 5, 6, 29, 26, 32, 13, 11, 30, 17, 31, 12] and many other softwares
available online). However, SVM cannot provide the confidence of class prediction since it
only learns the decision boundaries between diﬀerent classes. To cope with this problem,

Superfast-Trainable Multi-Class Probabilistic Classifier2

several post-processing methods have been developed for approximately computing theclass-posterior probability [25, 34].
On the other hand, logistic regression (LR) is a classiﬁcation algorithm that can naturally give the conﬁdence of class prediction since it directly learns the class-posterior
probabilities [15]. Recently, various eﬃcient algorithms for training LR models specialized in sparse data have been developed [22, 10].
Applying the kernel trick to LR as done in SVM, one can easily obtain a non-linear
classiﬁer with probabilistic outputs, called kernel logistic regression (KLR). Since the kernel matrix is often dense (e.g., Gaussian kernels), the state-of-the-art LR algorithms forsparse data are not applicable to KLR. Thus, in order to train KLR classiﬁers, standard
non-linear optimization techniques such as Newton’s method (more speciﬁcally, iteratively
reweighted least-squares) and quasi-Newton methods (for example, the Broyden-FletcherGoldfarb-Shanno (BFGS) method) seem to be commonly used in practice [15, 23]. Although the performance of these general-purpose non-linear optimization techniques has
been improved together with the evolution of computer environment in the last decade,computing the KLR solution is still challenging when the number of training samples islarge. The purpose of this paper is to propose an alternative probabilistic classiﬁcation
method that can be trained very eﬃciently.
Our proposed method is called the Least-Squares Probabilistic Classifier (LSPC). In
LSPC, we use a linear combination of Gaussian kernels centered at training points as amodel of class-posterior probabilities. Then we ﬁt this model to the true class-posterior
probability by least-squares1 . An advantage of this linear least-squares formulation isthat consistency is guaranteed without taking into account the normalization factor. In
contrast, normalization is essential in the maximum-likelihood LR formulation; otherwise
the likelihood tends to inﬁnity. Thanks to the simpliﬁcation brought by excluding thenormalization factor from the optimization criterion, we can compute the globally optimal
solution of LSPC analytically just by solving a system of linear equations.
Furthermore, we show that the use of a linear combination of kernel functions in
LSPC allows us to learn the parameters in a class-wise manner. This highly contributes tofurther reducing the computational cost particularly in multi-class classiﬁcation scenarios.
Through experiments, we show that LSPC is computationally much more eﬃcient than
KLR with comparable accuracy.

2

Least-squares Approach to Probabilistic Classification

In this section, we formulate the problem of probabilistic classiﬁcation and give a new
method in the least-squares framework.
1

A least-squares formulation has been employed for improving the computational eﬃciency of SVMs
[29, 26, 13]. However, these approaches deal with deterministic classiﬁcation, not probabilistic classiﬁcation.

Superfast-Trainable Multi-Class Probabilistic Classifier

2.1

3

Problem Formulation

Let X (⊂ Rd ) be the input domain, where d is the dimensionality of the input domain.
Let Y = {1, . . . , c} be the set of labels, where c is the number of classes. Let us consider
a joint probability distribution on X × Y with joint probability density p(x, y). Supposethat we are given n independent and identically distributed (i.i.d.) paired samples of
input x and output y:{(xi , yi ) ∈ X × Y}ni=1 .
The goal is to estimate the class-posterior probability p(y|x) from the samples{(xi , yi )}ni=1 . The class-posterior probability allows us to classify test sample x to class yb
with conﬁdence p(by |x):
yb := argmax p(y|x).
yLet us denote the marginal density of x by p(x) and we assume that it is strictly
positive:p(x) > 0 for all x ∈ X .
Then, by deﬁnition, the class-posterior probability p(y|x) can be expressed asp(y|x) =p(x, y)
.
p(x)

(1)

This expression will be utilized in the derivation of the proposed method below.

2.2

Linear Least-squares Fitting of Class-posterior Probability

Here we introduce our least-squares ﬁtting idea. We begin with the formulation forlearning the class-posterior probability p(y|x) as a function of both x and y, i.e., theclass-posterior probabilities for all classes are learned simultaneously. Then in Section 2.3,
we show that this simultaneous learning problem can be decomposed into independentclass-wise learning problems, which highly contributes to reducing the computational cost.
We model the class-posterior probability p(y|x) by the following linear model:
q(y|x; α) :=b∑

αℓ ϕℓ (x, y) = α⊤ ϕ(x, y),ℓ=1

where ⊤ denotes the transpose of a matrix or a vector,
α = (α1 , . . . , αb )⊤are parameters to be learned from samples, andϕ(x, y) = (ϕ1 (x, y), . . . , ϕb (x, y))⊤are basis functions such thatϕ(x, y) ≥ 0b for all (x, y) ∈ X × Y.

Superfast-Trainable Multi-Class Probabilistic Classifier

4

0b denotes the b-dimensional vector with all zeros and the inequality for vectors is appliedin the element-wise manner. We explain how the basis functions ϕ(x, y) are practicallychosen in Section 2.3.
We determine the parameter α in the model q(y|x; α) so that the following squared
error J is minimized:c ∫1∑J(α) :=(q(y|x; α) − p(y|x))2 p(x)dx
2 y=1
c ∫
c ∫
∑
1∑
2
=
q(y|x; α) p(x)dx −q(y|x; α)p(x, y)dx + Const.
2 y=1
y=1
1
= α⊤ Hα − h⊤ α + Const.,
2
where we used Eq.(1). The b × b matrix H and the b-dimensional vector h are deﬁned as
H :=
h :=

c ∫
∑
y=1
c ∫
∑

ϕ(x, y)ϕ(x, y)⊤ p(x)dx,
ϕ(x, y)p(x, y)dx.

y=1

H and h contain the expectations over unknown densities p(x) and p(x, y), so we approximate the expectations by sample averages. Then we have
c
n
1 ∑∑
c
H :=
ϕ(xi , y)ϕ(xi , y)⊤ ,n y=1 i=1
n
1∑
b
h :=
ϕ(xi , yi ).
n i=1

Now our optimization criterion is formulated as
[
]
⊤
1 ⊤c
⊤
b
b := argmin α Hα − h α + λα α ,α2α∈Rbwhere a regularizer λα⊤ α (λ > 0) is included for regularization purposes. Taking thederivative of the above objective function and equating it to zero, we see that the solutionb can be obtained just by solving the following system of linear equations.
αbc + λI b )α = h,(H(2)

b is given anawhere I b denotes the b-dimensional identity matrix. Thus, the solution αlytically asbc + λI b )−1 h.
b = (HαSuperfast-Trainable Multi-Class Probabilistic Classifier

5b is a conditional probability, we round up
In order to assure that the solution q(y|x; α)negative outputs to zero [35] and renormalize the solution. Consequently, our ﬁnal solution
is expressed asb ⊤ ϕ(x, y)max(0, αpb(y|x) = ∑c,(3)b ⊤ ϕ(x, y ′ ))y ′ =1 max(0, αWe call the above method Least-Squares Probabilistic Classifier (LSPC). LSPC can be
regarded as an application of a density ratio estimation method called the unconstrained
Least-Squares Importance Fitting (uLSIF) [18] to probabilistic classiﬁcation. Thus all thetheoretical properties of uLSIF such as consistency, the rate of convergence, and numericalstability [19, 20] may be directly translated into the current context.

2.3

Basis Function Design

A naive choice of basis functions ϕ(x, y) would be a kernel model, i.e., for some kernel
function K ′ ,nc ∑∑(y ′ )
αℓ K ′ (x, xℓ , y, y ′ ),(4)q(y|x; α) =y ′ =1 ℓ=1

which contains cn parameters. For this model, the computational complexity for solving
Eq.(2) is O(c3 n3 ).
Here we propose to separate input x and output y, and use the delta kernel for y (asin KLR):nc ∑∑(y ′ )αℓ K(x, xℓ )δy,y′ ,q(y|x; α) =y ′ =1 ℓ=1

where K is a kernel function for x and δy,y′ is the Kronecker delta:{
1 if y = y ′ ,δy,y′ =
0 otherwise.
This model choice actually allows us to speed up the computation of LSPC signiﬁcantly
since all the calculations can be carried out separately in a class-wise manner. Indeed,
the above model for class y is expressed asq(y|x; α) =

n
∑

(y)

αℓ K(x, xℓ ).

(5)

ℓ=1

c becomes block-diagonal, as illustrated in Figure 1(a). Thus we only
Then the matrix H
need to train a model with n parameters separately c times for each class y, by solving
the following equation:e (y) ,
c ′ + λI n )α(y) = h(HSuperfast-Trainable Multi-Class Probabilistic Classifier

Kernels
for class 1H

Kernels
for class 2

6

Samplesin class 1

Kernels
for class 3HSamplesin class 2

Samplesin class 3

(1)

′


H

′


H


H

(2)

′

(3)

H

(a) Model (5)(b) Model (7)c for model (5) and model (7). The number of classes is
Figure 1: Structure of matrix H
c = 3. Suppose training samples {(xi , yi )}ni=1 are sorted according to label y. Colored
blocks are non-zero and others are zeros. For model (5) consisting of c sets of n basis
c becomes block-diagonal (with common block matrix Hc ′ ), andfunctions, the matrix H
thus training can be carried out separately for each block. For model (7) consisting of csets of ny basis functions, the size of the target block is further reduced.
e (y) is the n-dimensional vector deﬁned asc ′ is the n × n matrix and hwhere Hn1∑′bHℓ,ℓ′ :=K(xi , xℓ )K(xi , xℓ′ ),
n i=1
n
1∑
(y)
e
K(xi , xℓ )δy,yi .
hℓ :=
n i=1

c ′ is common to all y, we only need to compute (Hc ′ + λI n )−1 once. Then the
Since H
computational complexity for obtaining the solution is O(n3 + cn2 ), which is smaller thanthe case with general kernel model (4). Thus this approach would be computationallyeﬃcient when the number of classes c is large.
Here, we further propose to reduce the number of kernels in model (5). To this end,
we focus on a kernel function K(x, x′ ) that is “localized”. Examples of such localized
kernels include the popular Gaussian kernel [28]:)(∥x − x′ ∥2′.
(6)
K(x, x ) = exp −2σ 2
Our idea is to reduce the number of kernels by locating the kernels only at samples
belonging to the target class:q(y|x; α) =

ny
∑
ℓ=1

(y)

(y)

αℓ K(x, xℓ ),

(7)

7

OutputSuperfast-Trainable Multi-Class Probabilistic Classifier

Input

Figure 2: Heuristic of reducing the number of basis functions—locate Gaussian kernels
only at the samples of the target class.
(y) nywhere ny is the number of training samples in class y, and {xi }i=1
is the training input
samples in class y.
The rationale behind this model simpliﬁcation is as follows (Figure 2). By deﬁnition,
the class-posterior probability p(y|x) takes large values in the regions where samples inclass y are dense; conversely, p(y|x) takes smaller values (i.e., close to zero) in the regions
where samples in class y are sparse. When a non-negative function is approximated by aGaussian kernel model, many kernels may be needed in the region where the output of
the target function is large; on the other hand, only a small number of kernels would beenough in the region where the output of the target function is close to zero. Following
this heuristic, many kernels are allocated in the region where p(y|x) takes large values,
which can be achieved by Eq.(7).
This model simpliﬁcation allows us to further reduce the computational cost since thec is further reduced, as illustrated in Figure 1(b). In
size of the target blocks in matrix H
order to learn the ny -dimensional parameter vector

α(y) = (α1 , . . . , αn(y)
)⊤
y
(y)

for each class y, we only need to solve the following system of ny linear equations:c(Hcwhere H

(y)

(y)

b ,
+ λI ny )α(y) = h
(y)

(8)

b (y) is the ny -dimensional vector deﬁned asis the ny × ny matrix and h
n
1∑
(y)
(y)
(y)
b
Hℓ,ℓ′ :=K(xi , xℓ )K(xi , xℓ′ ),
n i=1

(9)

ny
1∑
(y)
(y)
(y)
b
K(xi , xℓ ).
hℓ :=
n i=1

b (y) be the solution of Eq.(8). Then our ﬁnal solution is given byLet α∑ny (y)(y)max(0, ℓ=1αbℓ K(x, xℓ ))pb(y|x) = ∑c.
∑ny′ (y′ )(y ′ )bℓ K(x, xℓ ))ℓ=1 αy ′ =1 max(0,(10)

Superfast-Trainable Multi-Class Probabilistic Classifier

Input:

Output:8

Labeled training samples {(xi , yi )}ni=1(y) ny(equivalently, {xi }i=1
for class y = 1, . . . , c),Gaussian width σ, and regularization parameter λ;Class-posterior probability pb(y|x);

for y = 1, . . . , c

(
)
n
(y) 2
(y) 2
∑
∥
+
∥x
−
x
∥
1
∥x
−
x
′
i
i
ℓ
ℓ
b (y)′ ←−exp −Hfor ℓ, ℓ′ = 1, . . . , ny ;
ℓ,ℓ
n i=1
2σ 2
(
)
ny
(y)
(y) 2
∑
∥x
−
x
∥
1
(y)
ℓ
b
hℓ ←−
exp − ifor ℓ = 1, . . . , ny ;
n i=1
2σ 2
b (y) and obtain αc (y) + λI ny )α(y) = hb (y) ;
Solve linear equation (Hend(
max 0,

pb(y|x) ←−

c
∑

(

ny
∑

(y)

α
bℓ expℓ=1
ny ′max 0,

y ′ =1(∑(y)∥x − xℓ ∥2
−
2σ 2

(

(y ′ )
αbℓ expℓ=1))(y ′ )∥x − xℓ ∥2
−
2σ 2

)) ;

Figure 3: Pseudo code of LSPC for simpliﬁed model (7) with Gaussian kernel (6).
For the simpliﬁed model (7), the computational complexity for obtaining the solutionis O(cn2y n)—when ny = n/c for all y, this is equal to O(c−1 n3 ). Thus this approach iscomputationally highly eﬃcient for multi-class problems.
A pseudo code of the simplest LSPC implementation for Gaussian kernels is summaR
rized in Figure 3. Its MATLAB⃝
implementation is available fromhttp://sugiyama-www.cs.titech.ac.jp/∼sugi/software/LSPC/

3

Experiments

In this section, we experimentally compare the performance of the following classiﬁcation
methods:• LSPC: LSPC with model (7).
• LSPC(full): LSPC with model (5).
• KLR: ℓ2 -penalized kernel logistic regression with Gaussian kernels. We used aR
MATLAB⃝
implementation included in the ‘minFunc’ package [27], which uses
limited-memory BFGS updates with Shanno-Phua scaling in computing the step
direction and a bracketing line-search for a point satisfying the strong Wolfe conditions to compute the step size.

Superfast-Trainable Multi-Class Probabilistic Classifier

9When we fed data to learning algorithms, the input samples were normalized in theelement-wise manner so that each element has mean zero and unit variance. The Gaussian
width σ and the regularization parameter λ for all the methods are chosen based on 2-foldcross-validation from1
m, 15 m, 21 m, 23 m, m, 32 m, 2m, 5m, 10m},
σ ∈ { 10λ ∈ {10−2 , 10−1.5 , 10−1 , 10−0.5 , 100 },where
m := median({∥xi − xj ∥}ni,j=1 ).

3.1

Illustrative Examples

First, we illustrate the behavior of each method using a toy dataset.
We set the dimension of the input space to d = 2 and the number of classes to c = 3.
We independently drew samples in each class from the following class-conditional sample
densities (see Figure 4):( [ ] [])
−2
1 0
p(x|y = 1) = N x;
,
,
0
0 1
])
( [ ] [
1 0
2,,
p(x|y = 2) = N x;
0 1
0
])
])
( [ ] [
( [ ] [
1
1
4 0
1 0
0
0
p(x|y = 3) = N x;
+ N x;
,
,
,
0 1
0 4
−3
2
2
2
where N (x; µ, Σ) denotes the Gaussian density with mean µ and covariance matrix Σ.
We set the class-prior probabilities p(y) as{
1/4 if y = 1, 2,
p(y) =
1/2 if y = 3,
and we set the number of training samples to n = 200. Generated samples are plotted in
Figure 5.
The true class-posterior probabilities p(y|x) (∝ p(x|y)p(y)), their estimates obtainedby LSPC, LSPC(full), and KLR are depicted in Figure 6. The plots show that all themethods approximate the true class-posterior probabilities well in the training region(say, [−5, 5]2 ). However, the output outside the training region is substantially diﬀerent
in LSPC and KLR. This is induced by the diﬀerence of the models—a linear combinationof Gaussian kernels is used in LSPC, while its exponent is used in KLR. Outside thetraining region, there is no kernel, and thus a linear combination of Gaussian kernels
takes values close to zero (note that the values are not exactly zero since Gaussian tails
extended from training regions remain everywhere); then typically one of the classes takes
a value close to one, and the others tend to zero outside the training regions. On the otherSuperfast-Trainable Multi-Class Probabilistic Classifierp(x|y=1)p(x|y=2)

10

5

0

x

0

0.1
(2)

(2)

0.1
x

(2)

x

p(x|y=3)

10

5

0.05
−5

−5

0
(1)
x

5

10

−10
−10

0

10

0.04

5

0.03

0

0.02

−5

0.01

0.05

−5

−10
−10

10

−5

0
(1)
x

5

−10
−10

0

10

−5

0
(1)
x

5

10

Figure 4: Illustrative examples. Class-conditional sample densities p(x|y).

Optimal classification

Classification by LSPC

10

10

8

8

6

6

4

4

2

2

0

0

−2

−2

Tr1
Te1
Tr2
Te2
Tr3
Te3

−4
−6
−8
−10
−10

−5

0

5

−4
−6
−8
10

−10
−10

−5

Classification by LSPC(full)
10

10

8

8

6

6

4

4

2

2

0

0

−2

−2

−4

−4

−6

−6

−8

−8

−10
−10

−10
−10

−5

0

0

5

10

5

10

Classification by KLR

5

10

−5

0

Figure 5: Illustrative examples. Training samples are plotted with ﬁlled symbols. Unﬁlled
symbols denote the classiﬁcation results based on the true class-posterior probabilities andtheir estimates obtained by LSPC, LSPC(full), and KLR.

Superfast-Trainable Multi-Class Probabilistic Classifier

True p(y=1|x)

True p(y=2|x)10

True p(y=3|x)

10
0.8

5

10
0.8

5

0
(1)
x

5

−10
−10

10

p−hat(y=1|x) by LSPC

(2)

0.2

−5

0
(1)
x

5

0.8

5

−10
−10

(2)

−5

0
(1)
x

5

10

10

−10
−10

0

(2)

−5

0
(1)
x

5

10

−10
−10

0

0.2

−5

0
(1)
x

5

10

(2)

0

0.6
0

0.4
−5

−10
−10

0

0.8

5

x

(2)

(2)

0.4

0.6
x

10

0.8

−5

0.8

0.4

5

1

p−hat(y=3|x) by KLR

0.6

0
(1)
x

0

10

5

−5

10

0.6

0.2

0.8
5

0.2

5

0

p−hat(y=2|x) by KLR

−5

0
(1)
x

5

10

0

−5

0.4

p−hat(y=1|x) by KLR

x

10

0

−5

10

−10
−10

1

x

(2)

(2)

x

x
5

−10
−10

0.2

0.6

0.2

0
(1)
x

0

0.8

5

0.4

−5

0.4
−5

p−hat(y=3|x) by LSPC(full)

0.6

−10
−10

0.8
0.6

0.2

10

0.8

−5

1

0

p−hat(y=2|x) by LSPC(full)
1

0

10

0.4

p−hat(y=1|x) by LSPC(full)

5

5

5

x

(2)

x

(2)

x

0

−5

10

10

0
(1)
x

0.6

0.2

5

−5

10

0.8

5

0.4
−5

0.2

p−hat(y=3|x) by LSPC
1

0.6
0

0
(1)
x

−10
−10

10

10

−5

0.4
−5

p−hat(y=2|x) by LSPC

10

−10
−10

0.6
0

0.4
−5

0.2

−5

0.8

5

x

(2)

0

0.4
−5

1

0.6
x

x

(2)

0.6
0

−10
−10

11

0.2

−5

0
(1)
x

5

10

0.4
−5

−10
−10

0.2
−5

0
(1)
x

5

10

Figure 6: Illustrative examples. The plots show the true class-posterior probabilitiesp(y|x), their estimates by LSPC, LSPC(full), and KLR from top to bottom, and y = 1, 2, 3
from left to right.

Superfast-Trainable Multi-Class Probabilistic Classifier

12

hand, KLR outputs values close to one outside the training region since exp(0) = 1; thenthey are normalized and thus are reduced to 1/c.
The classiﬁcation results based on the true class-posterior probabilities and their estimates obtained by LSPC, LSPC(full), and KLR are plotted in Figure 5. This shows that
all the method gave reasonable classiﬁcation results.

3.2

Performance Comparison

Next, we evaluate the classiﬁcation accuracy and computation time of each method using
the following multi-class classiﬁcation datasets taken from the LIBSVM web page [5]:• mnist: Input dimensionality is 717 and the number of classes is 10.
• usps: Input dimensionality is 256 and the number of classes is 10.
• satimage: Input dimensionality is 36 and the number of classes is 6.
• letter: Input dimensionality is 16 and the number of classes is 26.
We investigated the classiﬁcation accuracy and computation time of LSPC,
LSPC(full), and KLR. For given n and c, we randomly chose ny = ⌊n/c⌋ training samples
from each class y, where ⌊t⌋ is the largest integer not greater than t. In the ﬁrst set ofexperiments, we ﬁxed the number of classes c to the original number shown above, andchanged the number of training samples as n = 100, 200, 500, 1000, 2000. In the second
set of experiments, we ﬁxed the number of training samples to n = 1000, and changed
the number of classes c—samples only in the ﬁrst c classes in the dataset are used. The
classiﬁcation accuracy is evaluated using 100 test samples randomly chosen from eachclass. The computation time is measured by the CPU computation time required fortraining each classiﬁer when the Gaussian width and the regularization parameter chosen
by cross-validation were used.
The experimental results are summarized in Figure 7 and Figure 8. The left column
in Figure 7 shows that when n is increased, the misclassiﬁcation error for all the methods
tends to decrease, and LSPC, LSPC(full), and KLR performed similarly well. The rightcolumn in Figure 7 shows that when n is increased, the computation time tends to grow for
all the methods. LSPC is faster than KLR by two orders of magnitude. The left column
in Figure 8 shows that when c is increased, the misclassiﬁcation error tends to increase
for all the methods, and LSPC, LSPC(full), and KLR behaved similarly well. The rightcolumn in Figure 8 shows that when c is increased, the computation time of KLR tends
to grow, while that of LSPC is kept constant or even it tends to slightly decrease. This
happened because the number of samples in each class decreases when c is increased, andthe computation time of LSPC is governed by the number of samples in each class, notby the total number of samples (see Section 2.3).
Overall, the computation of LSPC was shown to be faster than that of KLR by ordersof magnitude, while LSPC and KLR were shown to be comparable to each other interms of the classiﬁcation accuracy. LSPC and LSPC(full) were shown to possess similarSuperfast-Trainable Multi-Class Probabilistic Classifier

Misclassification rate (mnist, c=10)13Training time (mnist, c=10)

4

10

0.4
LSPC
0.35

LSPC(full)
KLR

0.3

2

10

0.25
0.2
0

10

0.15
0.1
0.05
100

−2

200

500
n

1000

2000

Misclassification rate (usps, c=10)

10

100

200

500
n

1000

2000

2000

5000

Training time (usps, c=10)

4

10
0.22
0.2

2

10

0.18
0.16
0.14

0

10

0.12
0.1
0.08
0.06
100

−2

200

500

1000

2000

5000

10

100

200

500

1000
n

n
Misclassification rate (satimage, c=6)Training time (satimage, c=6)

4

10

0.4
0.38
0.36

2

10

0.34
0.32
0

10

0.3
0.28
0.26
100

−2

200

500
n

1000

2000

10

100

Misclassification rate (letter, c=26)
4

0.8

200

500
n

1000

2000

Training time (letter, c=26)

10

0.7
0.6

2

10

0.5
0.4
0

10

0.3
0.2
0.1
100

−2

200

500
n

1000

10

100

200

500

1000

n

Figure 7: Misclassiﬁcation rate (in percent, left) and computation time (in second, right)
as functions of the number of training samples n. From top to bottom, the graphs
correspond to the ‘mnist’, ‘usps’, ‘satimage’, and ‘letter’ datasets.

Superfast-Trainable Multi-Class Probabilistic Classifier

Misclassification rate (mnist, n=1000)

14

Training time (mnist, n=1000)

4

10

0.12
0.1

3

10

0.08

2

10
0.06

1

10
0.04

0

LSPC
LSPC(full)
KLR

0.02
0
2

4

6
c

8

10

−1

10

Misclassification rate (usps, n=1000)

10

2

4

6
c

8

10

8

10

Training time (usps, n=1000)

3

10
0.09
0.08

2

10

0.07
0.06

1

10

0.05
0.04
0

10

0.03
0.02

−1

0.01
2

4

6
c

8

10

10

Misclassification rate (satimage, n=1000)

2

4

6
c

Training time (satimage, n=1000)

3

10

0.35
0.3

2

10

0.25
0.2

1

10
0.15
0.1

0

10

0.05
−1

0
2

3

4
c

5

6

Misclassification rate (letter, n=1000)

10

2

3

4
c

5

6

20

26

Training time (letter, n=1000)

4

10

0.25
0.2

2

10
0.15
0.1

0

10
0.05
0
2

−2

5

10

15
c

20

26

10

2

5

10

15
c

Figure 8: Misclassiﬁcation rate (in percent, left) and computation time (in second, right)
as functions of the number of classes c. From top to bottom, the graphs correspond tothe ‘mnist’, ‘usps’, ‘satimage’, and ‘letter’ datasets.

Superfast-Trainable Multi-Class Probabilistic Classifier

15classiﬁcation performance, and thus a computationally eﬃcient version, LSPC, would be
more preferable in practice.

4

Discussion and Conclusion

Recently, various eﬃcient algorithms for computing the solution of logistic regression havebeen developed for high-dimensional sparse data [22, 10]. However, for dense data, using
standard non-linear optimization techniques such as Newton’s method or quasi-Newton
methods seem to be a common choice [15, 23]. The performance of these general-purpose
non-linear optimizers has been improved in the last decade, but computing the solution oflogistic regression for a large number of dense training samples is still a challenge problem.
In this paper, we proposed a simple probabilistic classiﬁcation algorithm called LeastSquares Probabilistic Classifier (LSPC). LSPC employs a linear combination of Gaussian
kernels centered at training points for modeling the class-posterior probability and theparameters are learned by least-squares. Notable advantages of LSPC are that its solution
can be computed analytically just by solving a system of linear equations and training
can be carried out separately in a class-wise manner. In experiments, we showed that
LSPC is faster than kernel logistic regression (KLR) in computation time by two ordersof magnitude, with comparable accuracy.
The computational eﬃciency of LSPC was brought by the combination of appropriate
model choice and loss function. More speciﬁcally, KLR uses a log-linear combination ofkernel functions and its parameters are learned by regularized maximum likelihood. In
this log-linear maximum likelihood formulation, normalization of the model is essential
to avoid the likelihood diverging to inﬁnity. Thus the likelihood function tends to becomplicated and numerically solving the optimization problem may be unavoidable. On
the other hand, in LSPC, we chose a linear combination of Gaussian kernel functions
for modeling the class-posterior probability and its parameters are learned by regularized
least-squares. This combination allowed us to obtain the solution analytically. When
Newton’s method (more speciﬁcally, iteratively reweighted least-squares) is used for learning the KLR model, a system of linear equations needs to be solved in every iteration
until convergence [15]. On the other hand, LSPC requires to solve a system of linear
equations only once.
We chose to separate the kernel for inputs and outputs, and adopted the delta kernel
for outputs (see Eq.(5)). This allowed us to perform the training of LSPC in a class-wisemanner. We showed that this contributes to reducing the training time particularly inmulti-class classiﬁcation problems. We note that this model choice is essentially the sameas that of KLR2 .
We further proposed to reduce the number of kernels when “localized” kernels suchas the Gaussian kernel (6) is used. Through the experimental evaluation in Section 3, we
found that this heuristic model simpliﬁcation does not degrade the classiﬁcation accuracy,
2
The number of parameters in LSPC with model (5) is cn, while the number of parameters in KLR is(c − 1)n since the normalization (‘sum-to-one’) constraint is incorporated in the training phase.

Superfast-Trainable Multi-Class Probabilistic Classifier16

but reduces the computation time.
It is straightforward to show that solutions for all regularization parameter values (i.e.,
the regularization path, see [9, 14]) can be computed eﬃciently in LSPC. Let us consider
c (y) (see Eq.(9)):
the eigendecomposition of the matrix Hny∑c (y) =
Hγℓ ψ ℓ ψ ⊤
ℓ ,

ℓ=1
y
c
where {ψ ℓ }ℓ=1
are the eigenvectors of H(y)b can be expressed as
the solution α

n

b
α

(y)

c
= (H

(y)

(y)

n

y
associated with the eigenvalues {γℓ }ℓ=1
. Then,−1 b (y)+ λI ny ) h

=

ny
∑
b ⊤ψ
h

ℓ

γℓ + λℓ=1

ψℓ.

b ⊤ ψ ℓ )ψ ℓ is common to all λ, we can compute the solution αb (y) for all λ eﬃciently
Since (hc (y) once in advance. Although eigendecomposition ofby eigendecomposing the matrix H
c (y) may be computationally slightly more demanding than solving a system of linearHequations of the same size, this approach would be useful, e.g., when computing thesolutions for various values of λ in the cross-validation procedure.
When ny is large, we may further reduce the computational cost and memory spaceby using only a subset of kernels.
q(y|x; α) =

by
∑

(y)

(y)

αℓ K(x, cℓ ),

ℓ=1
n
∑
(y)
(y)
b (y)′ = 1
H
K(xi , cℓ )K(xi , cℓ′ ),
ℓ,ℓ
n i=1
ny
1∑
(y)
(y) (y)b
K(xi , cℓ ),
hℓ =
n i=1
(y) b

(y) n

y
y
where by is a constant chosen to be smaller than ny and {cℓ }ℓ=1
is a subset of {xℓ }ℓ=1.
This would be a useful heuristic when a huge number of samples are used for training.
Another option for reducing the computation time when the number of samples is verylarge would be the stochastic gradient descent method [1]. That is, starting from some
initial parameter value, gradient descent is carried out only for a randomly chosen single
sample in each iteration. Since our optimization problem is convex, convergence to theglobal solution is guaranteed (in a probabilistic sense) by stochastic gradient descent.
We focused on using the delta kernel for class labels (see Section 2.3). We expect thatdesigning appropriate kernel functions for class labels would be useful for improving theclassiﬁcation performance, e.g., in the context of multi-task learning [4, 2, 21]. We will
pursue this direction in our future work.

Superfast-Trainable Multi-Class Probabilistic Classifier

17

Acknowledgments
The author thanks Dr. Ryota Tomioka, Mr. Jaak Simm, and Dr. Hirotaka Hachiya for
their valuable comments. This work was supported by AOARD, SCAT, and the JST
PRESTO program.

References[1] S. Amari, “Theory of adaptive pattern classiﬁers,” IEEE Transactions on Electronic
Computers, vol.EC-16, no.3, pp.299–307, 1967.
[2] B. Bakker and T. Heskes, “Task clustering and gating for Bayesian multitask learning,” Journal of Machine Learning Research, vol.4, pp.83–99, 2003.
[3] C.M. Bishop, Pattern Recognition and Machine Learning, Springer, New York, NY,
USA, 2006.
[4] R. Caruana, L. Pratt, and S. Thrun, “Multitask learning,” Machine Learning, vol.28,
pp.41–75, 1997.
[5] C.C. Chang and C.J. Lin, “LIBSVM: A library for support vector machines,” tech.
rep., Department of Computer Science, National Taiwan University, 2001. http://www.csie.ntu.edu.tw/∼cjlin/libsvm/.
[6] R. Collobert and S. Bengio., “SVMTorch: Support vector machines for large-scale
regression problems,” Journal of Machine Learning Research, vol.1, pp.143–160, 2001.
[7] C. Cortes and V. Vapnik, “Support vector networks,” Machine Learning, vol.20,
pp.273–297, 1995.
[8] R.O. Duda, P.E. Hart, and D.G. Stor, Pattern Classiﬁcation, Wiley, New York, 2001.
[9] B. Efron, T. Hastie, R. Tibshirani, and I. Johnstone, “Least angle regression,” The
Annals of Statistics, vol.32, no.2, pp.407–499, 2004.
[10] R.E. Fan, K.W. Chang, C.J. Hsieh, X.R. Wang, and C.J. Lin, “LIBLINEAR: A
library for large linear classiﬁcation,” Journal of Machine Learning Research, vol.9,
pp.1871–1874, 2008.
[11] R.E. Fan, P.H. Chen, and C.J. Lin, “Working set selection using second order information for training SVM,” Journal of Machine Learning Research, vol.6, pp.1889–
1918, 2005.
[12] V. Franc and S. Sonnenburg, “Optimized cutting plane algorithm for large-scale risk
minimization,” Journal of Machine Learning Research, vol.10, pp.2157–2192, 2009.

Superfast-Trainable Multi-Class Probabilistic Classifier

18[13] G.M. Fung and O.L. Mangasarian, “Multicategory proximal support vector machineclassiﬁers,” Machine Learning, vol.59, no.1–2, pp.77–97, 2005.
[14] T. Hastie, S. Rosset, R. Tibshirani, and J. Zhu, “The entire regularization path forthe support vector machine,” Journal of Machine Learning Research, vol.5, pp.1391–
1415, 2004.
[15] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning:
Data Mining, Inference, and Prediction, Springer, New York, 2001.
[16] T. Joachims, “Making large-scale SVM learning practical,” in Advances in Kernel Methods—Support Vector Learning, ed. B. Schölkopf, C.J.C. Burges, and A.J.
Smola, pp.169–184, The MIT Press, Cambridge, MA, 1999.
[17] T. Joachims, “Training linear SVMs in linear time,” ACM SIGKDD International
Conference on Knowledge Discovery and Data Mining (KDD2006), pp.217–226, 2006.
[18] T. Kanamori, S. Hido, and M. Sugiyama, “A least-squares approach to direct importance estimation,” Journal of Machine Learning Research, vol.10, pp.1391–1445,
Jul. 2009.
[19] T. Kanamori, T. Suzuki, and M. Sugiyama, “Condition number analysis of kernelbased density ratio estimation,” tech. rep., arXiv, 2009.
[20] T. Kanamori, T. Suzuki, and M. Sugiyama, “Theoretical analysis of density ratio
estimation,” IEICE Transactions on Fundamentals of Electronics, Communications
and Computer Sciences, vol.E93-A, no.4, pp.787–798, 2010.
[21] T. Kato, H. Kashima, M. Sugiyama, and K. Asai, “Conic programming for multitask learning,” IEEE Transactions on Knowledge and Data Engineering, vol.22, no.7,
pp.957–968, 2010.
[22] K. Koh, S.J. Kim, and S.P. Boyd, “An interior-point method for large-scale
l1 -regularized logistic regression,” Journal of Machine Learning Research, vol.8,
pp.1519–1555, 2007.
[23] T.P. Minka, “A comparison of numerical optimizers for logistic regression,” tech.
rep., Microsoft Research, 2007.
[24] J. Platt, “Fast training of support vector machines using sequential minimal optimization,” in Advances in Kernel Methods—Support Vector Learning, ed. B. Schölkopf,
C.J.C. Burges, and A.J. Smola, pp.169–184, The MIT Press, Cambridge, MA, 1999.
[25] J. Platt, “Probabilities for SV machines,” in Advances in Large Margin Classiﬁers,
ed. A.J. Smola, P.L. Bartlett, B. Schölkopf, and D. Schuurmans, The MIT Press,
Cambridge, MA, 2000.

Superfast-Trainable Multi-Class Probabilistic Classifier

19[26] R. Rifkin, G. Yeo, and T. Poggio, “Regularized least-squares classiﬁcation,” Advances in Learning Theory: Methods, Models and Applications, ed. J.A.K. Suykens,
G. Horvath, S. Basu, C. Micchelli, and J. Vandewalle, NATO Science Series III:
Computer & Systems Sciences, vol.190, Amsterdam, the Netherlands, pp.131–154,
IOS Press, 2003.
[27] M. Schmidt, minFunc, 2005. http://people.cs.ubc.ca/∼schmidtm/Software/
minFunc.html.
[28] B. Schölkopf and A.J. Smola, Learning with Kernels, MIT Press, Cambridge, MA,
2002.
[29] J.A.K. Suykens, T.V. Gestel, J.D. Brabanter, B.D. Moor, and J. Vandewalle, Least
Squares Support Vector Machines, World Scientiﬁc Pub. Co., Singapore, 2002.
[30] Y. Tang and H.H. Zhang, “Multiclass proximal support vector machines,” Journalof Computational and Graphical Statistics, vol.15, no.2, pp.339–355, 2006.
[31] C.H. Teo, Q. Le, A. Smola, and S.V.N. Vishwanathan, “A scalable modular convex
solver for regularized risk minimization,” ACM SIGKDD International Conferenceon Knowledge Discovery and Data Mining (KDD2007), pp.727–736, 2007.
[32] I. Tsang, J. Kwok, and P.M. Cheung, “Core vector machines: Fast SVM training
on very large data sets,” Journal of Machine Learning Research, vol.6, pp.363–392,
2005.
[33] V.N. Vapnik, Statistical Learning Theory, Wiley, New York, NY, USA, 1998.
[34] T.F. Wu, C.J. Lin, and R.C. Weng, “Probability estimates for multi-class classiﬁcation by pairwise coupling,” Journal of Machine Learning Research, vol.5, pp.975–
1005, 2004.
[35] M. Yamada, M. Sugiyama, G. Wichern, and J. Simm, “Improving the accuracyof least-squares probabilistic classiﬁers,” IEICE Transactions on Information and
Systems, 2011. submitted.