Distributed Bayesian Piecewise Sparse Linear Models

arXiv:1711.02368v1 [cs.AI] 7 Nov 2017

Masato Asahara
NEC System Platform Research Laboratories
masahara@nec-labs.com

Abstract—The importance of interpretability of machine
learning models has been increasing due to emerging enterprise
predictive analytics, threat of data privacy, accountability of
artificial intelligence in society, and so on. Piecewise linear
models have been actively studied to achieve both accuracy
and interpretability. They often produce competitive accuracy
against state-of-the-art non-linear methods. In addition, their
representations (i.e., rule-based segmentation plus sparse linear
formula) are often preferred by domain experts. A disadvantage of such models, however, is high computational cost for
simultaneous determinations of the number of “pieces” and
cardinality of each linear predictor, which has restricted their
applicability to middle-scale data sets. This paper proposes
a distributed factorized asymptotic Bayesian (FAB) inference
of learning piece-wise sparse linear models on distributed
memory architectures. The distributed FAB inference solves
the simultaneous model selection issue without communicating
O(N ) data where N is the number of training samples and
achieves linear scale-out against the number of CPU cores.
Experimental results demonstrate that the distributed FAB
inference achieves high prediction accuracy and performance
scalability with both synthetic and benchmark data.

I. I NTRODUCTION
The importance of interpretability and transparency of
machine learning models has been increasing due to emerging enterprise predictive analytics, threat of data privacy,
accountability of artificial intelligence in society, and so on.
In data mining and machine learning academic community,
the workshop named FAT/ML1 (fairness, accountability and
transparency in machine learning) has been held every year
since 2014. This momentum has grown by incorporating
legal aspects of machine learning agents (e.g., Symposium
on Machine Learning and the Law in NIPS20162 ). From
government point of view, European Union enforces GDPR3
(general data protection regulation) which requires that the
consequences of profiling (i.e., how models profile individuals) should be informed to the data subject. On the other
hand, interpretable models restrict model representations,
and the balance between interpretability and accuracy has
been important research topics for decades [1].
Piecewise linear models have been actively studied to
achieve both accuracy and interpretability, which include
from classical tree [2] or linear [3] ones to more advanced
1 http://www.fatml.org/
2 http://www.mlandthelaw.org/
3 http://www.eugdpr.org/

Ryohei Fujimaki
NEC Data Science Research Laboratories
rfujimaki@nec-labs.com

ones [4], [5], [6], [7]. They often produce competitive
accuracy against state-of-the-art non-linear methods on realworld datasets. In addition, their representations (i.e., rulebased segmentation plus sparse linear formula) are often
preferred by domain experts. To the best use of these
models, i.e., simple and accurate, we have to simultaneously
determine the number of “pieces” and cardinality of each
predictor. However, such simultaneous model selection is
essentially much more computationally demanding which
has restricted their applicability to middle-scale data sets.
Meanwhile, for analyzing very large scale data (a.k.a. Big
Data), the size of a feature matrix easily exceeds memory
capacity in a single computation node. Recent trends in
distributed computational platforms for large-scale machine
learning have been shifting from those based on distributed
file systems (e.g., Hadoop [8]) to those based on distributed
memories (e.g., Spark [9] and Parameter server [10]). While
Hadoop incurs substantial overhead due to the load of intermediate data from disks between computations, distributed
memory architectures are able to avoid the need for disk
access by storing data in memory across the computations.
Notably, Spark appears to be one of the most promising
platforms for enterprise data analytics, and many distributed
machine learning algorithms for Spark have recently been
developed [11], [12], [13], [14], [15], [16], [17], [18], [19],
[20].
This paper proposes a novel distributed algorithm for
learning piecewise linear models on distributed memory
architectures and an efficient implementation on Spark. Our
contributions are summarized as below.
Distributed Learning Algorithm: This paper develops a
distributed learning algorithm of piecewise linear models
with model selection. Our technical contributions are mainly
two-fold. First, our algorithm linearly scales over the number
of distributed workers and automates the model selection
problem by taking advantages of recently-developed two
techniques: 1) factorized asymptotic Bayesian hierarchical
mixture of experts (FAB/HME) [21] for model selection
of piecewise linear models and 2) median selection subset
aggregation estimator (MESSAGE) [22] for communicationefficient distributed feature selection. Second, the MESSAGE algorithm independently processes data in each
worker for communication efficiency and we observe this
yields a bias in factorized information criterion (FIC) [23],

[24], [25], [26], [27], the model selection criterion of
FAB/HME. We derive an asymptotic correction term of this
bias in FIC, which leads better feature selection of individual
local models.
Practical Design on Spark: This paper presents a design of
our algorithm on Spark that helps to fully utilize distributed
computation resources. A resilient distributed dataset (RDD)
is designed to perform iterative model optimization without shuffling data. Further, we show that the sample-wise
parallelization of our algorithm uses CPU resources much
more efficiently than expert-wise one. Our experimental
results demonstrate that our algorithm and design achieves
both high prediction accuracy and high scalability with both
synthetic and benchmark data compared to state-of-the-art
Spark machine learning libraries.
II. R ELATED W ORK
Piecewise linear models have actively been studied to
achieve both interpretability and accuracy. Such models include from classical ones such as decision trees [2], [28] and
Lasso [3] to more advanced models such as hierarchical mixture of experts (HME) [4], Bayesian treed linear models [29],
local supervised learning through space partitioning [5],
informative projection ensembles [7], supersparse linear
integer models [30] and so on. Optimization of piecewise
linear models is usually non-convex due to simultaneous
optimizations of partitions and local models. Partition-wise
linear models [6] addressed this issue by formulating it
as a structured-sparsity problem. FAB/HMEs [21] induce
sparseness both on tree structures and cardinalities of local
models, and fully automate simultaneous model selection
for learning piecewise linear models via FAB inference [23],
[26]. Jialei et al. [31] extended FAB/HMEs and incorporated
non-linearity in local predictors to gain better accuracy by
keeping a certain level of interpretability. Ribeiro et al. [32]
proposed to locally approximate non-linear models by linear
models for model agnostic interpretability. As far as we
know, the applications of sophisticated piecewise linear
models have been limited to middle scale datasets due to
their high computational costs.
Meanwhile, Spark [9] appears to be one of the most
promising platforms for distributed machine learning algorithms. There are a large quantity of researches to realize
distributed machine learning algorithms on Spark such as logistic regression [15], SVM [15], K-Means [16], LDA [17],
ADMM [18], dominant cluster detection [14], graph algorithms [13] and so on. Because of its high scale computing power, automation of hyper-parameter search on Spark
such as [33] is also an active research field. Furthermore,
Spark has gotten a lot of attention as a platform of deep
learning recently [19], [20]. These research outcomes are
continuously integrated with Spark as its machine learning
library called MLlib [12], [11]. As growing the proposals
of cutting-edge technologies, application field of machine

②

①

Figure 1: Example of Bernoulli Gating Node.

learning on Spark is spreading to industrial area including
electric power [34], telecommunication [35] and drug discovery [36]. In this movement, the importance of Spark as a
data science platform for KDD community is also growing
as we can find several tutorials held in the last several KDD
conferences [37], [38], [39]. Despite increasing attentions of
Spark-based high-scale machine learning, to the best of our
knowledge, there are little studies on distributed learning of
piecewise linear models (on Spark).
III. P RELIMINARY
Probabilistic piecewise sparse linear models, or
FAB/HME models [21], partition feature spaces using
gating functions and assign a sparse expert in each
partition. The FAB/HME model employs the Bernoulli
gating function as follows:
g(x, βi ) := gi U (ti − x[γi ]) + (1 − gi )U (x[γi ] − ti ), (1)
where βi = (gi , ti , γi ), gi ∈ [0, 1], U is a step function,
γi is the index w.r.t. the element of x ∈ RD , where D is
feature dimensionality, and ti ∈ R is an arbitrary value.
For example, when x[γi ] < ti , g(x, βi ) = gi and otherwise
g(x, βi ) = 1 − gi , as shown in Fig. 1.
Formally, the probabilistic models are defined as follows:
p(y|x, θ) =

E Y
X

ψg (x, i, j)p(y|x, φj ),

(2)

j=1 i∈Ej

where x ∈ RD is an observation variable, y is a target variable (y ∈ R for regression or y ∈ {1, −1} for classification),
θ = (β1 , . . . , βG ,
φ1 , . . . , φE ) represents model parameters, E and G are
the numbers of experts and gating functions, respectively.
Ej (j = 1, . . . , E) denotes the j-th expert index set and
contains all indices of the gating nodes on the unique path
(i,j)
from the root node to the j-th expert node. ψg (x) :=
ψ(g(x, βi ), i, j) is the probability on the i-th gate, and
Q
(i,j)
(x) expresses the probability of the j-th path,
i∈Ej ψg
where ψ(a, i, j) = a if the j-th expert is on the left subtree of the i-th gate, and 1 − a otherwise. The experts can
be either linear regression, i.e., p(y|x, φj ) = N (y|wjT x, σj2 )
where φj = (wj , σj2 ), or linear logistic regression, i.e.,
p(y|x, φj ) = 1/(1 + exp(−yφTj x)).
The latent variable related to the j-th path is defined by
ζj , where ζj = 1 if y|x is generated by the j-th node through

the j-th path, and ζj = 0 otherwise. The complete likelihood
is defined, then, as follows:
p(y N |ζ N , xN , φ) =

N Y
E
Y

(n)

p(y (n) |x(n) , φj )ζj ,

(3)

n=1 j=1

p(ζ N |xN , β) =

N Y
E Y
Y

(n)

ψg(i,j) (x(n) )ζj ,

(4)

n=1 j=1 i∈Ej

i∈Ej

The FAB/HME algorithm finds the best parameters and
models by maximizing the factorized information criterion
(FIC), which is an asymptotically accurate approximation
of Bayesian marginal log-likelihood [23], [26], derived as
follows:
G 
n h
X
F IC(y N , xN ) = max Eq log p(y N , ζ N |xN , θ)−
q

(t)

i=1

:::::

E
N X
N
 X
i
X
X
Dφj
Dβi
(n)
(n)
(n)
q(ζj ) −
`j q(ζj )
log
log
2
2
n=1
n=1
j=1
j∈Gi

::::::::::::::::::::::::::::::::::::::::::::::

o
+ Hq ,

(6)

where q is any distribution on ζ N , Hq is the entropy of
q, and Gi is the i-th gating index set and contains all
indices of the experts on the sub-tree of the i-th gating
(n)
(n)
node. `j is a scaling factor where `j
= 1/σj2 for
(n)
(n)
(n)
(n)
linear regression or `j = µj (1 − µj ) where µj =
1/(1 + exp(−y (n) φTj x(n) )) for logistic regression. This
optimization is conducted by alternating optimizations of q
and θ like expectation-maximization algorithm. Note that
Dφj is equivalent to the L0 -norm of φj , which induces
sparsity of the model. [21] has applied the forward-backward
greedy (FoBa) algorithm [40] for least square experts, which
has the tightest error bounds among state-of-the-art methods.
We can easily extend this idea to logistic regression by
applying the gradient FoBa algorithm [41], which has the
same theoretical error bounds for general smooth convex
functions.

(t)

Input: • Nβi , Nφj , θ(t)
Input: ? q (t) (ζ (n) ), θ(t)
Output: • F IC (t)

P
(t)
(n) (n) (n) (t)
1: ? LLw =
, ζ |x , θ )] −
n∈Nw Eq [log p(y

(t) (n)
(t) (n)
q (ζ ) log q (ζ )
for ∀n ∈ Nw and ∀w ∈ W
2:

where φ = (φ1 , . . . , φE ) and β = (β1 , . . . , βG ). Prediction
is executed by the j∗ -th expert that maximizes the gating
probability as follows:
Y
j∗ = arg max
ψg(i,j) (x).
(5)
j

Algorithm 1: Distributed FIC Computation

• F IC (t) =
PW
(t) PG
w=1 LLw −
i=1

(t)

(t)

Dβ
2

i

(t)

log Nβi −

PE

j=1

Dφ

j

2

like that on a shared memory architecture. The second challenge is to avoid huge communication overhead. Particularly
communications among worker nodes cause reallocation of
data (a.k.a. data shuffling in Spark) in distributed memory
architectures. Third, for the “map-reduce” computational
model which is one of the most popular modern distributed
memory computation models, balancing or equalizing CPU
loads is needed to minimize synchronization overhead.
Hereafter, we denote by D a dedicated server and by w ∈
W a worker node, where W is a set of worker nodes. We
assume that the training data are distributed on memories
in ∀w ∈ W 4 . The subscription w denotes the w-th worker
node. We denote by Nw a sample index set of the w-th
worker node. In Algorithm 1 ∼ Algorithm 5, ? and • are
executed on the worker nodes in parallel and the dedicated
server in serial, respectively.

A. Distributed FIC Computation
The distributed computation of FIC, which is necessary
for convergence determination, is described in Algorithm 1.
The FIC calculation consists of two parts: 1) the sum of
expected log-likelihood and 2) regularization. The former
requires xN , y N and q (t) (ζ N ), which are distributed on
the worker nodes, and hence each worker node computes
expected log-likelihoods of data in its memory, and then
only its summation needs to be collected by the dedicated
server, as shown in line 1 of Algorithm 1. This step requires
only communication of scalar values between the dedicated
server and individual worker nodes. The latter is computed
in the dedicated server, as shown in line 2 of Algorithm 1.

IV. D ISTRIBUTED FAB/HME A LGORITHM
This paper considers situations in which data size N is
much larger than data dimensionality D. The designation
of the distributed FAB/HME algorithm and architecture
involves the following three challenges. First, due to memory
capacity limitation, no single worker node can load entire
data, target, and variational distributions. This prohibits the
algorithm from employing a straightforward parallelization

(t)

log Nφj

B. Distributed E-step
The distributed computation of the E-step is described in
(n)
Algorithm 2. First, q(ζj ) for ∀n ∈ Nw is calculated on
4 The training data is stored in a distributed file system and loaded on
memory of workers in parallel.

Algorithm 2: Distributed E-step
(t−1)
(t−1)
Input: ? Nβi , Nφj , θ(t−1)
(t)
(t)
Output: • Nφj , Nβi
(n)
(t)
(t)
Output: ? q (t) (ζj ), Nφj , Nβi
(n,t)
1: ? Calculate qj
using (7) for ∀n ∈ Nw , ∀w ∈
P
P
(n,t)
2: ? Nφwj =
q
and Nβwi = j∈Gi Nφwj
n∈Nw j
P
P
(t)
(t)
w
w
3: • Nφ =
w∈W Nφj and Nβi =
w∈W Nβi
j

4:

Algorithm 3: Split Points Calculation
Input: • Tmax : number of the maximum split points
Output: ? Td for ∀d.
(n)
1: ? Calculate xmax
dw = maxn∈Nw {xd } and
(n)
xmin
dw = minn∈Nw {xd } for ∀d and ∀w.
2: • Calculate xmax
= maxw {xmax
d
dw } and
min
min
xd = maxw {xdw }.
max
3: • Calculate Td by splitting the range [xmin
] into
d , xd
Tmax bins with the equal width.

W.

• Eliminate irrelevant experts and gates using (8).

Algorithm 4: Distributed Bernoulli Gate Optimization
w ∈ W as follows:
(n,t)

qj

(n)

(t)

:= q (t) (ζj ) ∝

Y

(t−1)

ψg(i,j,t−1) (x(n) )p(y (n) |x(n) , φj

i∈Ej

exp

n X −D
βi
(t−1)

i∈Ej

+

2Nβi

(n,t−1) o

−Dφj `j

(t−1)

,

2Nφj

::::::::::::::::::::::::::::::

(7)
The expected numbers of samples on experts and gates,
denoted by Nφwj and Nβwi , are then collected by the dedicated server, and, therefore, only two scalar values are
(t)
(t)
communicated. In line 3, Nφj and Nβi are computed and
globally shared on the dedicated server. It is known that
exponentiated regularization (the waved part of (7)) derived
from FIC (the waved part of (6)) eliminates redundant latent
variables through EM iterations [23]. Such “non-effective”
experts are eliminated from the model as follows:
(n,t)

qj

(t)

(n,t)

= 0 if Nφj < δ otherwise qj

(t)

/Qj ,

(8)

(t)
and Qj are a threshold value and a
PE (n,t)
for j=1 qj
= 1. This shrinkage

where δ
normalization
constant
process automatically determines piecewise space partitioning structures
in a Bayesian-principled fashion.
C. Distributed M-step: Bernoulli Gates
The gate optimization requires a set of split candidate
points, i.e., Bd = {ti |discretized domain of xd }, for each
dimension d of x. To consistently aggregate distributed
calculations, all worker nodes must share Bd for ∀d. For
this purpose, Bd is computed at the beginning of the FABEM algorithm on the basis of Algorithm 3. First, on each
(n)
worker node, the maximum and minimum values of xd
are computed and then collected by the dedicated server.
Here a vector of the size D is transferred from each worker
node to the dedicated server. Bd is computed on the basis
of lines 2 and 3 in Algorithm 3. B = {Bd } is distributed
to all worker nodes. This process communicates a matrix of
the size DTmax .
The distributed gate optimization is summarized in Algorithm 4. First, some intermediate statistics are computed in

)

(n)

Input: ? Td , Nβi , q (t) (ζj )
(t)
Output: ? βi
1: for Tγi where γi = 1, . . . , D do
R
2:
? Calculate ρL
w (γi , ti ) and ρw (γi , ti ) for ∀w ∈ W
using (9) and (10).
3: end for
4: for Tγi where γi = 1, . . . , D do
P
5:
• Calculate ρL (γi , ti ) = w∈W ρL
w (γi , ti ) and
P
R
R (t) (t)
ρ (γi , ti ) = w∈W ρw (γi , ti ).
6:
• Calculate ξi (γi , ti ) where
(t)
gi (γi , ti ) = (ρL (γi , ti ) + ρR (γi , ti ))/Nβi .
7: end for
(t) (t)
8: • Calculate γi , ti = arg maxγi ,ti ξi (γi , ti ), and
(t)
(t) (t)
gi = gi (γi , ti ).

a distributed manner as follows (at line 1-3):
X
X (n,t)
ρL
qj
w (γi , ti ) =

(9)

n∈Nw ∩Tsi j∈GiL

ρR
w (γi , ti ) =

X

X

(n,t)

qj

(10)

n∈Nw ∩Tli j∈GiR

where Tli , Tsi are the sets of samples whose γi -th dimension
is larger or smaller than ti , GiL contains all indices of the
expert nodes on the left sub-tree of the i-th gating node,
and GiR is similarly defined for the right sub-tree of the iR
th gating node. The two matrices w.r.t. ρL
w and ρw , each of
which has DTmax elements, are collected by the dedicated
server. Then, on the dedicated server, the i-th gate parameter
(t)
βi is computed as described in lines 4-8 in Algorithm 4,
(t)
and βi , containing three scalar values, is distributed to all
worker nodes.
D. Distributed M-step: Sparse Experts
For optimizing experts in the M-step, we have to distribute L0 regularized optimization, to which well-studied
approaches using distributed gradient or proximal methods
are not applicable. We address this issue by applying a
recently-developed median selection subset aggregation estimator (MESSAGE) algorithm [22].

Algorithm 5: Distributed Sparse Experts Optimization
(n)

(t−1)

Input: ? q (t) (ζj ), φj
(t)
Output: ? φj
?
•
3: ?
4: •
1:
2:

(t+ 1 )

Calculate φjw 2 using (11) for ∀w ∈ W.
Perform the majority voting using (12).
(t)
Estimate φjw using (13).
(t)
(t)
Compute φj by aggregating φjw using (14).

N

|W|

log p(y (n) |x(n) , φj ) −

X (n,t−1) (n,t)
Dφj
log
`j
qj ,
2
n=1

it causes different bias in FIC by reusing the first term (i.e.,
the first term is evaluated |W| times on a single data set).
In order to avoid this bias, we consider an asymptotic
approximation of the first term as follows5 :
X (n,t)
qj (log p(y (n) |x(n) , φ∗j ) − log p(y (n) |x(n) , φ̂j )) ≈ −Dφj
n∈Nw

(11)
N
X

Dφj
(n,t−1) (n,t)
−
`j
qj .
log
2
n=1

The detailed derivation
of (11) is discussed in the next
(t+ 12 )
sub-section. φjw
is then once collected by the dedicated
server, and majority voting is performed to determine a
feature set as follows (line 2):
X (t+ 1 )
|W|
Fj = {d |
[φjw 2 ]d ≥
}.
(12)
2
w∈W

The feature set Fj is then distributed to all workers, and
parameter estimation is performed using only features in Fj
as follows (line 3):
X (n,t)
(t)
φjw = arg max
qj
log p(y (n) |x(n) , φj (Fj )), (13)
φj

(n,t)

qj

(16)

n∈Nw

− (|W| − 1)Dφj

X
n∈Nw

The distributed optimization of sparse experts is described
in Algorithm 5. First, feature selection is performed using
the FoBa algorithm on each worker node as follows (line 1):
X (n,t)
(t+ 1 )
qj
log p(y (n) |x(n) , φj )
φjw 2 = arg max |W|
φj

The first (loss) and the second (regularizer) terms are O(N )
and O(log N ), respectively.
PN Therefore,
P in distributed M-step,
if we simply replace
n=1 by
n∈Nw , it changes the
balance between the first and the second terms. On the other
hand, if we rescale the first term as follows:

n∈Nw

where φj (Fj ) means that parameters not included in Fj are
(t)
fixed to zero. Finally, φjw is again collected by the dedicated
server, and the weight aggregation is performed to estimate
the weight vector as follow (line 4):
1 X (t)
(t)
φj =
φjw ,
(14)
|W|
w∈W

(t)
φj

where
is distributed to all worker nodes. Note that the
majority voting and the weight aggregation are lightweight
computations in comparison to the L0 optimization processes, and their executions on the dedicated server do not
affect the parallel performance.
E. Correction of FIC in M-step
This subsection explains the derivation of (11). In (6), the
terms related to expert optimization can be summarized as
follows:
N
N
X
X
Dφj
(n,t−1) (n,t)
(n,t)
qj
log p(y (n) |x(n) , φj ) −
log
`j
qj .
2
n=1
n=1
(15)

N
X

(n,t)

qj

(log p(y (n) |x(n) , φ∗j ) − log p(y (n) |x(n) , φ̄j )) ≈ −Dφj

n=1

(17)
where φ∗ is the true parameter and φ̂j and φ̄j are
P
(n,t)
the maximizer of
log p(y (n) |x(n) , φj ) and
n∈Nw qj
PN
(n,t)
(n) (n)
log p(y |x , φj ), respectively. By taking
n=1 qj
into account that Nw can be considered as a random subset
of N , we have:
N
X

(n,t)

qj

log p(y (n) |x(n) , φj )

(18)

n=1

= |W|ENw [

X

(n,t)

qj

log p(y (n) |x(n) , φ∗j )]

n∈Nw

where ENw is the expectation over the randomness on Nw .
By combining (15), (17) and (18), we have (11).
V. E FFICIENT D ESIGN ON S PARK
This section describes an efficient design of the distributed
FAB/HME algorithm on Spark. Hereinafter we refer to the
distributed FAB/HME as dFAB and the original FAB/HME
as sFAB.
A. RDD Structure and Execution Flow
A resilience distributed dataset (RDD) is the base in
Spark distributed computation, on whose elements all distributed computations on Spark are performed, i.e., Spark
processes each element of an RDD in parallel. An RDD
is an immutable and partitioned collection of records and
can only be created from persistent data or other RDDs
via transformations. A standard RDD design might assign
one data instance to an element of the RDD. However,
major computations of FAB/HME algorithms rely on matrix
computations, and the RDD for dFAB has to be designed to
5 The derivation can be obtained in a similar manner of the derivation of
Akaike’s information criterion [42].

X1 Y1 Q1

L1

B

Le

B

RDD(t-1)

(t+1)-th
FIC calculation

FIC calculation

・
・
・
e-th element
w-th partition

Xe

Ye Qe

Distributed M-step
(Gate Opt.)
RDD(t-1)
w/ Le(t)

・
・
・
・
・
・

Distributed
E-step

Distributed M-step
(Experts Opt.)

RDD(t)

FoBa and
Majority Voting

(t-1)-th
model parameters
(t-1)

β

Figure 2: RDD structure of dFAB.

(t)

β

(t+1)-th E-step
^F

(t-1)

φ

FoBa and
Weight Average

Serialized in closures
and sent to executors

(t)

φ

Figure 3: Execution flow of dFAB.

i∈Ej

Fig. 2 illustrates the RDD structure in dFAB. One partition
contains multiple elements that are data units to which map
functions are applied, and one element of w-th partition
consists of a tuple of (Xe , Ye , Qe , Le , B), where e is the
index of an element and w corresponds to the worker node
index w7 . Fig. 3 illustrates the execution flow of dFAB. Let
us denote the RDD after the t-th EM iteration by RDD(t) ,
(t)
(t)
whose element is (Xe , Ye , Qe , Le , B).
• In the t-th FIC calculation, the driver process invokes
map transformation with a closure executing Step 1 of
Algorithm 1, which involves model parameters β (t−1)
and φ(t−1) as its parameters. Here, each executor com(t)
putes Le as an intermediate outcome, and each RDD
(t−1)
(t−1)
, B) to
, Le
element is updated from (Xe , Ye , Qe
(t)
(t−1)
, Le , B). Note that, if the eliminated ex(Xe , Ye , Qe
perts (Step 4 of Algorithm 1) are identified in the driver
(t−1)
s are updated in the executors by simply
process, Qe
setting the corresponding column to be zero in order to
avoid reallocation of the RDD itself.
• In the t-th E-step, the driver process invokes a map
transformation with a closure executing Steps 1 and 2
in the Algorithm 2. Note that we avoid recalculation
(t)
of expert-wise likelihood by storing Le as a part of
the RDD element. Here, each RDD element is updated
(t−1)
(t)
(t)
(t)
from (Xe , Ye , Qe
, Le , B) to (Xe , Ye , Qe , Le , B).
Through the FIC calculation and E-step, we obtain RDD(t)
without re-allocating or shuffling any large scale data.
• The t-th M-step is performed on the basis of Algorithm 4
and 5. This does not change RDD(t) .
6 Although we can use mapPartition function to collect a part of
training data from a partition of the RDD, mapPartition passes an
iterator of samples to the closure, which causes many iterator accesses and
significant computational overhead.
7 For the evaluations in this paper we set the number of elements so that
it equals the number of workers.

# active experts

10 4
Log-scaled frequency

10 3
10 2
10 1
10 0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

process them efficiently6 . For convenience, let us introduce
a few notations. Let Xe , Ye , and Qe be matrices whose
elements are x(n) for ∀n ∈ Ne , y (n) for ∀n ∈ Ne , and
q(ζ (n) ) for ∀n ∈ Ne , respectively, where Ne ⊂ Nw . Note
(t)
that each element of Qe is q (t) (ζ (n) ) for ∀n ∈ Ne . Also,
(t)
Le is a |Ne | × E matrix, and the (n, j)-th element of Le
is:
Y
(t−1)
ψg(i,j,t−1) (x(n) )p(y (n) |x(n) , φj
).
log

Time for optimizing an expert (sec)

35
30
25
20
15
10
50 10 20 30 40 50 60 70
FAB EM iterations

Figure 4: Histogram of Figure 5: Number of active
time for experts optimiza- experts over FAB EM itertion.
ations.

B. Advantage of Sample-wise Parallelization over Expertwise Parallelization
This subsection discusses advantages of sample-wise parallelization (i.e., Algorithm 5) over expert-wise parallelization. Although the later one might seem to be a natural
choice for parallelizing the FAB/HME learning algorithm,
since expert optimizations can be processed independently,
we observe two issues. First, as Fig. 4 shows, the distribution
of processing time for expert optimization is heavy tailed.
The optimization with long process time has the large
(t−1)
number of assigned instances, Nφj , which varies across
experts through the EM iterations. This implies (and we
have empirically observed) that the loads of CPU cores get
unbalanced, and fast ones have to wait for slow ones, as
shown in Fig. 6b, resulting in poor CPU usage. Second, The
number of CPU cores that expert-wise parallelization can
use is bounded by the number of experts. In the FAB/HME
learning procedure, irrelevant experts are automatically eliminated, and the number of experts decreases through the
iterations. Fig. 5 illustrates an example of decreasing numbers of active experts over FAB EM iterations on a simple
simulation dataset. In this experiment, we used a machine
with 16 CPU cores. The algorithm utilized only half CPU
cores, however, after the number of active experts became
8 at the 20-th step.
Fig. 6a illustrates our sample-wise parallelization of expert optimization. Our sample-wise parallelization addresses
both of the issues. For the unbalance of loads across the experts, the time complexity of expert optimization is approximately linear w.r.t. sample size, and hence the sample-wise

gate1

gate3

CPU time

E3

E2

E4

E4

E3

E3

E2

E2

E4
E3

CPU time

gate2
E4

E2

gate4

E3
E4

E2

expert1

E1

E1

E1

E1

E1

core1

core2

core3

core4

(a) Sample-wise

core1

core2

core3

expert2

expert5

core4

(b) Expert-wise

Figure 6: Two options for parallelizing expert optimization
in FAB M-step.

parallelization can equally distribute computational loads.
In addition, for degree of parallelization, since each task
processes optimization for all experts, shrinkage decrease in
the number of experts accelerates overall computation rather
than decreasing CPU usage.

expert3

expert4

Figure 7: The true tree structure of the artificial data.

intermediately-generated elements of an RDD and to result
in a need for their re-computation. To deal with the issues,
the dFAB makes the checkpoints at every 20 iterations,
which number has been determined on the basis of empirical
knowledge from our experiments.
VI. E XPERIMENTS AND D ISCUSSIONS

C. Practical Implementation Tips
All the distributed operations of dFAB are implemented
as chains of RDD transformations. The dFAB first converts
the training data into an RDD and repeatedly transforms one
RDD into another RDD to calculate the variational distribution q and the model parameters θ. A driver process invokes
the ? processes of dFAB on executors by means of such
RDD transformers as map and reduce, and the calculation
results are collected in the driver process via collect or
count functions. Since all the distributed operations are
RDD transformations, dFAB is able to obtain the benefits
of Spark, e.g., sophisticated task/job scheduling [43] and
RDD fault tolerance.
A driver process and executors exchange the model parameters θ via dynamically-generated closures. A driver
distributes θ by generating a closure passed to a transformation function of Spark. The dFAB implements the logic
of dFAB as first-class parameterized functions. By passing
θ as a parameter of the parameterized function, the function
generates another closure which calculates its logic with
parameter θ. The dFAB passes the closure, derived from
θ to transformation functions, that executes dFAB logic on
executors.
Checkpointing of intermediate RDDs is important w.r.t.
the correct execution and the performance improvement.
dFAB is a kind of EM algorithm, which sometimes requires
long iterations. The dFAB executes the EM iterations by
chaining RDD operations. A long chain of RDD operations
increases the stack size of remote procedure calls and induces stack overflow at worst because Spark has to maintain
the lineage of the operations for keeping fault tolerance.
Furthermore, it is likely to cause garbage collection for

We used Spark 1.6.0, and 8 worker nodes ran Spark executors on YARN [44]. Each server and worker node employed
two Intel Xeon E5-2640 v3 processors (16 physical cores),
a 256 GB memory, and a 1 TB 7.2K RPM HDD, and they
were connected via a 1 Gbps Ethernet. Observed and target
variables were standardized in advance.
A. Comparison with the Serial Algorithm
We first demonstrate that dFAB improves its execution
speed with no negative impact on accuracy by comparing
it with the serial algorithm, which is abbreviated as sFAB,
using an artificial regression data set used in the original
FAB/HME paper [21]. The true tree structure illustrated in
Fig. 7 has 5 experts, and each expert uses 10–20 features.
On the i-th gating node, gi was fixed to 1, γi and ti
were randomly selected from [1, D] and [0, 1], respectively.
On the j-th expert, the non-zero elements of wj were
randomly sampled from [0, 1]. xN and y N were sampled
from Uniform[0, 1] and N (y (n) |wjT x(n) , 0.1), where N =
1, 000, 000, 5, 000, 000, 10, 000, 000 and D = 100. We
employed δ = 5 × 10−9 F IC (t−1) (termination condition)
and  = 3 × 10−2 N (shrinkage threshold). The number of
initial experts was 32 (5-depth symmetric tree).
Fig. 8 illustrates the prediction accuracy of dFAB and
sFAB over |W|. We observed:
• The RMSE values of dFAB were smaller than those
of sFAB. This might have been caused by the samplewise parallelization which reduced the training variance
as the original paper of the MESSAGE algorithm has also
reported [22].
• The RMSE of dFAB slightly increased with an increasing
|W| with relatively small N . This might be because the

0.30

N = 10 6

N = 5 × 10 6

N = 10 7

We compared the predictive accuracy of dFAB on 3 regression (gas sensor array (CO), gas sensor array (methane)
and household power consumption) and 2 classification data
sets (HIGGS and HEPMASS) in UCI repository [46], with
other distributed machine learning algorithms implemented
in Spark MLlib. The brief task abstracts of the data sets are
described here;

0.25
0.20
RMSE

B. Benchmark Evaluation

0.15
0.10

•

|W|

sFAB

1
16
32
64
128

1
16
32
64
128

0.00

1
16
32
64
128

0.05

dFAB

•

Figure 8: Accuracy comparison between dFAB and sFAB.

number of samples on each executor becomes insufficient.
However, the difference in accuracy became negligible
with larger N , which was over 5,000,000 in this experimental condition.
Fig. 9 shows the processing times of dFAB and sFAB over
|W| (i.e., number of CPU cores). We implemented parallel
gate optimization and expert-wise parallelization in sFAB
by OpenMP [45]. From these results, we observed that:
•

•

•

While EM iteration time for sFAB with 16 CPU cores
was as much as 73.8% less than that with 1 CPU core,
that for dFAB with 16 CPU cores was as much as 90.4%
less than the same8 .
dFAB with 128 CPU cores (i.e., |W| = 128) showed
an average time per EM iteration and a total processing
time as much as 98.4% and 98.5% less than those for
sFAB with 1 CPU core, respectively. In general, larger N
achieves a better improvement of time per EM iteration
with large |W |. A larger N leads to a longer CPU processing time for each worker and, as a result, diminishes
overhead in Spark execution, such as initial setup overhead
for executors and RDDs and I/O wait for intermediate
RDD checkpointing.
There existed a limitation in performance improvement
over |W| depending on the value of N . When N =
1, 000, 000, the performance improvement for EM iteration where |W| = 128 was worse than that where
|W| = 64. On the other hand, when N = 5, 000, 000
and 10, 000, 000, the performance improvement for EM
iteration where |W| = 128 was greater than that where
|W| = 64. This implies we should determine an adequate
number of partitions, i.e., |W|, on the basis of the size of
training data.

8 Due to the memory limitation, sFAB where |W| = 16 could not finish
where N = 10, 000, 000.

•

•

The task of gas sensor array (CO) and (methane) data sets
is multivariate regression with multiple responses based
on chemical sensing system which consists of an array
of 16 metal-oxide gas sensors and an external mechanical
ventilator to simulate the biological respiration cycle.
The task of household power consumption data set is to
predict the amount of power consumption every minutes
of a house, which contains three meter’s measurements
of electric power consumption in one household with a
one-minute sampling rate over a period of almost 4 years.
The task of HIGGS data set is a classification problem
to distinguish between a signal process which produces
Higgs bosons and a background process which does not,
which contains 42 features in total. We used 27 features
of all due to the limitation of the dFAB implementation.
The task of HEPMASS data set is to separate particleproducing collisions from a background source. This data
set produced by Monte Carlo simulations contains 28
features.

From Spark MLlib we chose ElasticNet [47], decision tree
(DecTree), random forests (RF) for both classification and
regression as baseline algorithms. Note that RF is less
interpretable and we evaluated it as a state-of-the-art distributed non-linear model. We used 2-loop cross validation
with 10-fold outer loops for evaluating test prediction error
and 3-fold inner loops for parameter selection. Note that
FAB/HMEs do not need 2-loop cross validation, so dFAB
does not execute the inner loop. We employed δ = 10−4
(termination condition) and  = N × 10−2 (shrinkage
threshold). The number of initial experts was 8 (3-depth
symmetric tree).
Table I summarizes test RMSEs for regression data sets
and classification errors for classification data sets. For all
three regression data sets, dFAB outperformed the others
while RF performed better than dFAB on the HIGGS data
for classification. For the HIGGS data set, dFAB generated
FAB/HMEs with 2-5 active experts whose cardinalities were
5-14, that were much more interpretable than the models
learned by RF, which consists of 300 trees with 5-depth. In
summary, the results indicate that 1) dFAB achieved better
predictive accuracy than that of other distributed algorithms
implemented in Spark MLlib, and 2) it achieved competitive
accuracy with more interpretable models than non-linear
models of Spark MLlib.
Next, we compared the reduction in processing time for

0.6
0.4

|W|

sFAB

1
16
32
64
128

0.0

1
16
32
64
128

0.2

dFAB

N = 5 × 10 6

N = 10 7

1
16
32
64
128

0.8

N = 10 6

1
16
32
64
128

N = 10 7

1
16
32
64
128

N = 5 × 10 6

Normalized time per FAB iteration

N = 10 6

1
16
32
64
128

Normalized processing time

1.0

1.0
0.8
0.6
0.4
0.2
0.0

sFAB

|W|

(a) Processing time until convergence.

dFAB

(b) Time per FAB EM iteration.

Figure 9: Processing times of dFAB and sFAB. The Y-axis represents normalized time by that of sFAB with 1 CPU core.
Table I: Test prediction errors (RMSEs for regression and 0-1 errors for classification). The numbers shown in parentheses
are standard deviations. The best and second best methods are highlighted in bold and bold italic faces, respectively.
data
gas sensor array (CO)
gas sensor array (methane)
household power consumption
HIGGS
HEPMASS

sample
4,208,261
4,178,504
2,075,259
11,000,000
7,000,000

dim
16
16
3
27
28

dFAB with that for other algorithm implementations of
Spark MLlib over |W| as shown in Fig. 10. We observed
that:
• dFAB outperformed other implementations with gas sensor array (methane), HIGGS and HEPMASS.
• With gas sensor array (CO), while dFAB outperformed
ElasticNet and DecTree, its improvement did not reach
the level of RF. The difference from gas sensor array
(methane) data set was caused by the difference of convergence time derived from the difference of training data,
even if the scale of data was similar to each other.
• dFAB had a limitation to scale out its computing performance on the household power consumption data. This is
because the data had only three features, which resulted
in running the FAB EM optimization processes quite
fast, and the overhead of Spark computation became a
bottleneck.
In summary, these results reveal the performance scalability
of dFAB against state-of-the-art machine learning implementation on Spark.
Finally, we compared the processing time for dFAB with
that for other algorithm implementations of Spark MLlib.
In these experiments we set |W| = 128 for all algorithm
implementations. Note that all processing time except dFAB

dFAB
0.542(0.01)
0.548(0.00)
0.524(0.01)
0.335(0.01)
0.156(0.00)

ElasticNet
0.597(0.00)
0.600(0.00)
0.531(0.00)
0.358(0.00)
0.163(0.00)

DecTree
0.587(0.01)
0.561(0.00)
0.529(0.00)
0.337(0.00)
0.167(0.00)

RF
0.576(0.00)
0.553(0.00)
0.655(0.00)
0.317(0.00)
0.175(0.00)

includes 3-fold inner loop. As shown in Table II, dFAB
averagely takes longer time than other algorithms because
dFAB solves non-convex optimization problems to achieve
both high interpretability and accuracy. However, as shown
in Fig. 10, the execution performance scalability of dFAB
is better than others. This property of dFAB alleviates the
performance disadvantage of dFAB because we can easily
reduce the execution time of dFAB by just adding more
worker nodes, which is not an unthinkable solution in this
cloud computing era.
C. Detailed Analysis on Resource Utilization
This section demonstrates that dFAB can empirically
utilize computing resources efficiently. Fig. 11 illustrates the
resource consumption of dFAB which runs with the artificial
data set and where N = 10, 000, 000. In the results, we
observed:
• CPU utilization of workers was roughly 60% during the
execution. Some decrease in CPU utilization occurred due
to the checkpointing RDDs. However, this did not severely
degrade the performance of the execution.
• The CPU core of the dedicated server for the driver
process was mostly idle, and the process was not a
bottleneck in the execution.

dFAB

dFAB

dFAB

ElasticNet

ElasticNet

ElasticNet

DecTree

DecTree

DecTree

RF
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Normalized processing time

RF
0.0

|W| = 16

|W| = 16

|W| = 64

|W| = 128

(a) gas sensor array (CO).

|W| = 64

|W| = 128

|W| = 16

(b) gas sensor array (methane).

dFAB

dFAB

ElasticNet

ElasticNet

DecTree

DecTree

RF
0.0

RF
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Normalized processing time

0.2 0.4 0.6 0.8 1.0
Normalized processing time

RF
0.0

0.2 0.4 0.6 0.8 1.0
Normalized processing time

|W| = 16

|W| = 64

|W| = 128

|W| = 128

(c) household power consumption.

0.2 0.4 0.6 0.8 1.0
Normalized processing time

|W| = 16

(d) HIGGS.

|W| = 64

|W| = 64

|W| = 128

(e) HEPMASS.

Figure 10: Comparison of the reduction of processing time between dFAB and other algorithm implementations of Spark
MLlib. Processing time is normalized by time of the same algorithm where |W| = 16
Table II: Processing time (seconds) with |W| = 128.
data
gas sensor array (CO)
gas sensor array (methane)
household power consumption
HIGGS
HEPMASS

•
•

•

sample
4,208,261
4,178,504
2,075,259
11,000,000
7,000,000

There were no spikes in memory utilization that implies
dFAB managed the memory utilization well.
Some spikes in disk and network access occurred due to
the checkpointing RDDs on HDFS. On the other hand, the
average of disk and network accesses during the iterations
without the checkpointing were less than several hundreds
of KB/s.
The driver process did not issue notable disk accesses.
This means that the execution of the driver was fully completed on memory and was not a performance bottleneck.
Further, the driver process issued network traffic at the rate
of tens of MB/s at the start of the execution. This network
transmission arose from the distribution process for the
training data. Once the process finished, the driver process
transferred insignificant amounts of data with workers,
which is less than hundreds of KB/s. We can diminish the
network transfer for training data by simple modification
of the implementation to make the workers read the data
from HDFS directly.

In summary, our results demonstrate that there are no serious
performance bottlenecks in dFAB and that efficient resource

dim
16
16
3
27
28

dFAB
1,145
379
289
1,396
1,307

ElasticNet
1,845
1,827
49
196
117

DecTree
55
68
27
75
74

RF
361
475
206
968
973

utilization resulted in good performance scalability in dFAB.
D. Scale-out Scalability Evaluation
To evaluate the scale-out (or horizontal scaling) scalability
of dFAB we conducted experiments on a cluster with tens to
hundreds virtual instances in Amazon Web Services (AWS).
We prepared two scaled-out clusters which have a large
number of worker nodes whose CPU has a small number
of cores; one consists of c3.large instances (denoted as
c3.large cluster) and another one consists of r3.large
instances (denoted as r3.large cluster). For comparison,
we also prepared two scaled-up clusters which have a
small number of worker nodes whose CPU has a large
number of cores; one consists of c3.8xlarge instances
(denoted as c3.8xlarge cluster) and another one consists of r3.8xlarge instances (denoted as r3.8xlarge
cluster). The hardware specification of them is described in
Table III. This evaluation used gas sensor array (methane)
data set.
Fig. 12a and 12b show average processing time and time
per EM iteration of dFAB on the clusters, respectively. The
experimental results show that the execution speed of dFAB

CPU usage (%)

worker3

100
80
60
40
20
0
100
80
60
40
20
0

worker6

0

10000

worker1

100
80
60
40
20
0

worker4

100
80
60
40
20
0
100
80
60
40
20
0

worker7

0

10000

worker2

100
80
60
40
20
0

worker5

100
80
60
40
20
0
100
80
60
40
20
0

Elapsed time (sec)

Active memory (GB)

driver

100
80
60
40
20
0

worker8

0

10000

driver

120
100
80
60
40
20
0

worker3

120
100
80
60
40
20
0
120
100
80
60
40
20
0

Disk I/O (MB/s)

worker3

60
50
40
30
20
10
0
60
50
40
30
20
10
0

worker6

0

10000

worker1

60
50
40
30
20
10
0

worker4

60
50
40
30
20
10
0
60
50
40
30
20
10
0

worker7

0

Elapsed time (sec)

10000

worker2

60
50
40
30
20
10
0

worker5

60
50
40
30
20
10
0
60
50
40
30
20
10
0

Read

worker7

0

10000

worker2

120
100
80
60
40
20
0

worker5

120
100
80
60
40
20
0
120
100
80
60
40
20
0

Elapsed time (sec)

worker8

0

10000

(b) Memory utilization.

Network I/O (MB/s)

driver

60
50
40
30
20
10
0

120
100
80
60
40
20
0

10000

(a) CPU utilization.

worker4

120
100
80
60
40
20
0

worker6

0

worker1

120
100
80
60
40
20
0

worker8

0

10000

Write

(c) Transfer rate of disk I/Os.

driver

50
40
30
20
10
0

worker3

50
40
30
20
10
0
50
40
30
20
10
0

50
40
30
20
10
0

10000

0

10000

worker5

50
40
30
20
10
0

worker7

Elapsed time (sec)

worker2

50
40
30
20
10
0

worker4

50
40
30
20
10
0

worker6

0

worker1

50
40
30
20
10
0

50
40
30
20
10
0

Read

worker8

0

10000

Write

(d) Transfer rate of network I/Os.

Figure 11: Resource consumption of dFAB.
Table III: Hardware specification of AWS instances.
Type
c3.large
c3.8xlarge
r3.large
r3.8xlarge

CPU cores
2
32
2
32

Memory (GiB)
3.75
60
15.25
244

Storage (GB)
2x16 SSD
2x320 SSD
1x32 SSD
2x320 SSD

on c3.large and r3.large clusters (scaled-out clusters)
is around 2 times faster than that of c3.8xlarge and
r3.8xlarge clusters (scaled-up clusters)9 . dFAB achieved
good scale-out scalability on c3.large and r3.large
clusters despite scaled-out architectures tend to increase the
network communications due to huge number of worker
nodes. This scale-out scalability is obtained by the design of
dFAB which reduces the network traffic as much as possible
to avoid performance degradation.
In summary, our experimental results demonstrate that the
design of dFAB algorithm and implementation possesses the
9 Even with same total number of CPU cores, the execution time
on c3.large and r3.large clusters was shorter than that on
c3.8xlarge and r3.8xlarge clusters. This is because good scaled-out
scalability of dFAB exploits the performance advantages of the hardware
architecture such as an automatically CPU clock-up technology and the
large total bandwidth of memory access.

Networking Perf.
Moderate
10 Gb
Moderate
10 Gb

Intel
Intel
Intel
Intel

Processor
Xeon E5-2680
Xeon E5-2680
Xeon E5-2670
Xeon E5-2670

v2
v2
v2
v2

Clock (GHz)
2.8
2.8
2.5
2.5

capability of high performance scalability on a scaled-out
cluster.
VII. C ONCLUSIONS
In this paper we have proposed the distributed FAB/HME
algorithm scalable algorithm to learn interpretable and accurate piecewise sparse linear models from Big Data. By taking
advantages of FAB inference and the MESSAGE algorithm
with a novel asymptotic bias correction of FIC, dFAB realizes fully automated model selection with linearly scale-out
capability over the data size. Further, we presented a design
of dFAB on Spark, the rising distributed memory computation platform. Our RDD for dFAB enables us to fully utilize
CPU without needing to shuffle data during optimization
processes. Our experimental results have demonstrated that
dFAB achieves high prediction accuracy and performance

Processing time (sec)

Time per FAB iteration (sec)

30

1000

r3.8xlarge
r3.large
c3.8xlarge
c3.large

800
600
400
200
00

50

100 150
# CPU cores

200

250

(a) Processing time until convergence.

r3.8xlarge
r3.large
c3.8xlarge
c3.large

25
20
15
10
5
00

50

100
150
# CPU cores

200

250

(b) Time per FAB EM iteration.

Figure 12: Processing time comparison between scaled-out clusters (c3.large and r3.large) and scaled-up clusters
c3.8xlarge and r3.8xlarge).

scalability for both synthetic and public data. One of important future work is to explore the possibility of rebuilding
the dFAB algorithm so as to make it asynchronous, in order
to improve the efficiency of computing resource utilization
and accelerate algorithm execution.
R EFERENCES
[1] A. A. Freitas, “Comprehensible Classification Models: A
Position Paper,” SIGKDD Explor. Newsl., vol. 15, no. 1, pp.
1–10, 2014.

[11] X. Meng, J. Bradley, B. Yavuz, E. Sparks, S. Venkataraman, D. Liu, J. Freeman, D. Tsai, M. Amde, S. Owen,
D. Xin, R. Xin, M. J. Franklin, R. Zadeh, M. Zaharia, and
A. Talwalkar, “MLlib: Machine Learning in Apache Spark,”
J. Mach. Learn. Res., vol. 17, no. 1, pp. 1235–1241, 2016.
[12] E. R. Sparks, A. Talwalkar, V. Smith, J. Kottalam, X. Pan,
J. Gonzalez, M. J. Franklin, M. I. Jordan, and T. Kraska,
“MLI: An API for Distributed Machine Learning,” in IEEE
ICDM, 2013, pp. 1187–1192.

[2] L. Breiman, J. Friedman, R. Olshen, and C. Stone, Classification and Regression Trees. Wadsworth and Brooks, 1984.

[13] J. E. Gonzalez, R. S. Xin, A. Dave, D. Crankshaw, M. J.
Franklin, and I. Stoica, “GraphX: Graph Processing in a
Distributed Dataflow Framework,” in USENIX OSDI, 2014,
pp. 599–613.

[3] R. Tibshirani, “Regression Shrinkage and Selection via the
Lasso,” J. Royal Statist. Soc., Series B, vol. 58, no. 1, pp.
267–288, 1996.

[14] L. Chu, S. Wang, S. Liu, Q. Huang, and J. Pei, “ALID:
Scalable Dominant Cluster Detection,” in VLDB, 2015.

[4] M. I. Jordan and R. A. Jacobs, “Hierarchical Mixtures of
Experts and the EM Algorithm,” Neural Comput., vol. 6,
no. 2, pp. 181–214, 1994.

[15] C. Y. Lin, C. H. Tsai, C. P. Lee, and C. J. Lin, “Large-scale
logistic regression and linear support vector machines using
spark,” in IEEE Big Data, 2014, pp. 519–528.

[5] J. Wang and V. Saligrama, “Local Supervised Learning
through Space Partitioning,” in NIPS, 2012, pp. 91–99.

[16] S. Gopalani and R. Arora, “Comparing Apache Spark and
Map Reduce with Performance Analysis using K-Means,”
IJCA, vol. 113, no. 1, pp. 8–11, 2015.

[6] H. Oiwa and R. Fujimaki, “Partition-wise Linear Models,” in
NIPS, 2014, pp. 3527–3535.
[7] M. Fiterau and A. Dubrawski, “Projection Retrieval for Classification,” in NIPS, 2012, pp. 3023–3031.
[8] Apache Hadoop. [Online]. Available: http://hadoop.apache.
org/
[9] M. Zaharia, M. Chowdhury, T. Das, A. Dave, J. Ma, M. McCauly, M. J. Franklin, S. Shenker, and I. Stoica, “Resilient
Distributed Datasets: A Fault-Tolerant Abstraction for InMemory Cluster Computing,” in USENIX NSDI, 2012, pp.
15–28.
[10] M. Li, D. G. Andersen, J. W. Park, A. J. Smola, A. Ahmed,
V. Josifovski, J. Long, E. J. Shekita, and B.-Y. Su, “Scaling
Distributed Machine Learning with the Parameter Server,” in
USENIX OSDI, 2014, pp. 583–598.

[17] Z. Qiu, B. Wu, B. Wang, C. Shi, and L. Yu, “Collapsed
Gibbs Sampling for Latent Dirichlet Allocation on Spark,”
in BIGMINE, 2014, pp. 17–28.
[18] S. Dhar, C. Yi, N. Ramakrishnan, and M. Shah, “ADMM
based scalable machine learning on Spark,” in IEEE Big Data,
2015, pp. 1174–1182.
[19] P. Moritz, R. Nishihara, I. Stoica, and M. I. Jordan,
“SparkNet: Training deep networks in Spark,” in ICLR, 2016.
[20] H. Kim, J. Park, J. Jang, and S. Yoon, “DeepSpark: SparkBased Deep Learning Supporting Asynchronous Updates and
Caffe Compatibility,” in NIPS, 2016.
[21] R. Eto, R. Fujimaki, S. Morinaga, and H. Tamano, “FullyAutomatic Bayesian Piecewise Sparse Linear Models,” in
AISTATS, 2014, pp. 238–246.

[22] X. Wang, P. Peng, and D. B. Dunson, “Median Selection
Subset Aggregation for Parallel Inference,” in NIPS, 2014,
pp. 2195–2203.
[23] R. Fujimaki and S. Morinaga, “Factorized asymptotic
bayesian inference for mixture modeling,” in AISTATS, 2012,
pp. 400–408.
[24] R. Fujimaki and K. Hayashi, “Factorized Asymptotic
Bayesian Hidden Markov Models,” in ICML, 2012, pp. 799–
806.
[25] K. Hayashi and R. Fujimaki, “Factorized Asymptotic
Bayesian Inference for Latent Feature Models,” in NIPS,
2013, pp. 1214–1222.
[26] K. Hayashi, S. Maeda, and R. Fujimaki, “Rebuilding
Factorized Information Criterion: Asymptotically Accurate
Marginal Likelihood,” in ICML, 2015, pp. 1358–1366.
[27] C. Liu, L. Feng, R. Fujimaki, and Y. Muraoka, “Scalable
Model Selection for Large-Scale Factorial Relational Models,” in ICML, vol. 37, 2015, pp. 1227–1235.
[28] S. K. Murthy, S. Kasif, and S. Salzberg, “A System for
Induction of Oblique Decision Trees,” Journal of Artificial
Intelligence Research, vol. 2, pp. 1–32, 1994.
[29] H. A. Chipman, E. I. George, and R. E. Mcculloch, “Bayesian
Treed Generalized Linear Models,” Bayesian Statistics, vol. 7,
2003.

[38] M. Armbrust, D. Bateman, R. Xin, and M. Zaharia, “Introduction to Spark 2.0,” KDD Tutorial, 2016.
[39] J.-M. Agosta, D. GuhaThakurta, R. H. M. Inchiosa, S. Kumar,
V. Paunic, H. Zhang, and M. Zhao, “Scalable R on Spark,”
KDD Tutorial, 2016.
[40] T. Zhang, “Adaptive Forward-Backward Greedy Algorithm
for Learning Sparse Representations,” IEEE Trans. on Information Theory, vol. 57, no. 7, pp. 4689–4708, 2011.
[41] J. Liu, R. Fujimaki, and J. Ye, “Forward-Backward Greedy
Algorithms for General Convex Smooth Functions over A
Cardinality Constraint,” in ICML, 2014, pp. 503–511.
[42] H. Akaike, “A new look at the statistical model identification,”
IEEE Trans. on Automatic Control, vol. 19, no. 6, pp. 716–
723, 1974.
[43] M. Zaharia, D. Borthakur, J. Sen Sarma, K. Elmeleegy,
S. Shenker, and I. Stoica, “Delay Scheduling: A Simple
Technique for Achieving Locality and Fairness in Cluster
Scheduling,” in Proc. of ACM SIGOPS European Conf.
on Comp. Sys., 2010, pp. 265–278. [Online]. Available:
http://doi.acm.org/10.1145/1755913.1755940
[44] V. K. Vavilapalli, A. C. Murthy, C. Douglas, S. Agarwal,
M. Konar, R. Evans, T. Graves, J. Lowe, H. Shah, S. Seth,
B. Saha, C. Curino, O. O’Malley, S. Radia, B. Reed, and
E. Baldeschwieler, “Apache Hadoop YARN: Yet Another
Resource Negotiator,” in ACM SoCC, 2013.

[30] B. Ustun, S. Tracà, and C. Rudin, “Supersparse Linear Integer
Models for Predictive Scoring Systems,” in AAAI Conference
on Late-Breaking Developments in the Field of Artificial
Intelligence, 2013, pp. 128–130.

[45] OpenMP. [Online]. Available: http://openmp.org/

[31] J. Wang, R. Fujimaki, and Y. Motohashi, “Trading Interpretability for Accuracy: Oblique Treed Sparse Additive
Models,” in ACM KDD, 2015, pp. 1245–1254.

[47] H. Zou and T. Hastie, “Regularization and variable selection
via the Elastic Net,” Journal of the Royal Statistical Society,
Series B, vol. 67, pp. 301–320, 2005.

[32] M. T. Ribeiro, S. Singh, and C. Guestrin, “”Why Should I
Trust You?”: Explaining the Predictions of Any Classifier,”
in ACM KDD, 2016, pp. 1135–1144.
[33] E. R. Sparks, A. Talwalkar, D. Haas, M. J. Franklin, M. I.
Jordan, and T. Kraska, “Automating Model Search for Large
Scale Machine Learning,” in ACM SoCC, 2015, pp. 368–380.
[34] J. Zheng and A. Dagnino, “An initial study of predictive
machine learning analytics on large volumes of historical data
for power system applications,” in IEEE Big Data, 2014, pp.
952–959.
[35] M. A. Alsheikh, D. Niyato, S. Lin, H. p. Tan, and Z. Han,
“Mobile big data analytics using deep learning and apache
spark,” IEEE Network, vol. 30, no. 3, pp. 22–29, 2016.
[36] D. Harnie, A. E. Vapirev, J. K. Wegner, A. Gedich, M. Steijaert, R. Wuyts, and W. D. Meuter, “Scaling Machine Learning for Target Prediction in Drug Discovery using Apache
Spark,” in Proc. of IEEE/ACM CCGrid, 2015, pp. 871–879.
[37] J. Shanahan and L. Dai, “Large Scale Distributed Data
Science using Apache Spark,” KDD Tutorial, 2015.

[46] A. Frank and A. Asuncin, “UCI machine learning repository,”
http://archive.ics.uci.edu/ml/, 2010.

