Identifying Consistent Statements about
Numerical Data with Dispersion-Corrected
Subgroup Discovery

arXiv:1701.07696v2 [cs.AI] 23 Apr 2017

Mario Boley∗†
mboley@mmci.uni-saarland.de
Bryan R. Goldsmith†
goldsmith@fhi-berlin.mpg.de
Luca M. Ghiringhelli†
ghiringhelli@fhi-berlin.mpg.de

Jilles Vreeken∗
jilles@mpi-inf.mpg.de

Abstract
Existing algorithms for subgroup discovery with numerical targets do not optimize the error or target variable dispersion of the groups they find. This often
leads to unreliable or inconsistent statements about the data, rendering practical
applications, especially in scientific domains, futile. Therefore, we here extend the
optimistic estimator framework for optimal subgroup discovery to a new class of
objective functions: we show how tight estimators can be computed efficiently for
all functions that are determined by subgroup size (non-decreasing dependence),
the subgroup median value, and a dispersion measure around the median (nonincreasing dependence). In the important special case when dispersion is measured
using the mean absolute deviation from the median, this novel approach yields a
linear time algorithm. Empirical evaluation on a wide range of datasets shows that,
when used within branch-and-bound search, this approach is highly efficient and
indeed discovers subgroups with much smaller errors.

1. Introduction
Subgroup discovery is a well-established KDD technique (Klösgen (1996); Friedman and
Fisher (1999); Bay and Pazzani (2001); see Atzmueller (2015) for a recent survey) with
applications, e.g., in Medicine (Schmidt et al, 2010), Social Science (Grosskreutz et al,
∗
†

Max Plank Institute for Informatics and Saarland University, Saarbrücken, Germany
Fritz Haber Institute of the Max Planck Society, Berlin

1

(a) optimizing coverage times median shift;
σ0 (·) ≡ (a(·) ≥ 6) ∧ (c2 (·) > 0) ∧ (c6 (·) <
v. high); subgroup median 0.22, subgroup
error 0.081

(b) optimizing
dispersion-corrected
variant;
σ1 (·) ≡ (a(·) ∈ [8, 12])∧(c2 (·) > low)∧(c6 (·) <
v. high) ∧ (r(·) > v. low); subgroup median 0.23, subgroup error 0.028

Figure 1. To gain an understanding of the contribution of long-range van der Waals interactions
(y-axis; above) to the total energy (x-axis; above) of gas-phase gold nanoclusters, subgroup
discovery is used to analyze a dataset of such clusters simulated ab initio by density functional
theory (Goldsmith et al, 2017); available features describe nanocluster geometry and contain,
e.g., number of atoms a, fraction of atoms with i bonds ci , and radius of gyration r. Here, similar
to other scientific scenarios, a subgroup constitutes a useful piece of knowledge if it conveys a
statement about a remarkable amount of van der Waals energy (captured by the group’s central
tendency) with high consistency (captured by the group’s dispersion/error); optimal selector σ0
with standard objective has high error and contains a large fraction of gold nanoclusters with a
target value below the global median (0.13) (a); this is not the case for selector σ1 discovered
through dispersion-corrected objective (b), which therefore can be more consistently stated to
describe gold nanoclusters with high van der Waals energy.

2010), and Materials Science (Goldsmith et al, 2017). In contrast to global modeling,
which is concerned with the complete characterization of some variable defined for a
given population, subgroup discovery aims to detect intuitive descriptions or selectors
of subpopulations in which, locally, the target variable takes on a useful distribution.
In scientific domains, like the ones mentioned above, such local patterns are typically
considered useful if they are not too specific (in terms of subpopulation size) and indicate
insightful facts about the underlying physical process that governs the target variable.
Such facts could for instance be: ‘patients of specific demographics experience a low
response to some treatment’ or ‘materials with specific atomic composition exhibit a
high thermal conductivity’. For numeric (metric) variables, subgroups need to satisfy
two criteria to truthfully represent such statements: the local distribution of the target
variable must have a shifted central tendency (effect), and group members must be
described well by that shift (consistency). The second requirement is captured by the
group’s dispersion, which determines the average error of associating group members
with the central tendency value (see also Song et al, 2016).
Despite all three parameters—size, central tendency, and dispersion—being impor-

2

tant, the only known approach for the efficient discovery of globally optimal subgroups,
branch-and-bound search (Webb, 1995; Wrobel, 1997), is restricted to objective functions that only take into account size and central tendency. That is, if we denote by Q
some subpopulation of our global population P then the objective functions f currently
available to branch-and-bound can be written as
f (Q) = g(|Q|, c(Q))

(1)

where c is some measure of central tendency (usually mean or median) and g is a function
that is monotonically increasing in the subpopulation size |Q|. A problem with all such
functions is that they inherently favor larger groups with scattered target values over
smaller more focused groups with the same central tendency. That is, they favor the
discovery of inconsistent statements over consistent ones—surprisingly often identifying
groups with a local error that is almost as high or even higher than the global error
(see Fig. 1 for an illustration of this problem that abounded from the authors’ research
in Materials Science). Although dispersion-corrected objective functions that counterbalance size by dispersion have been proposed (e.g., ‘t-score’ by Klösgen, 2002 or ‘mmad’
by Pieters et al, 2010), it remained unclear how to employ such functions outside of
heuristic optimization frameworks such as greedy beam search (Lavrač et al, 2004) or
selector sampling (Boley et al, 2012; Li and Zaki, 2016). Despite often finding interesting
groups, such frameworks do not guarantee the detection of optimal results, which can
not only be problematic for missing important discoveries but also because they therefore
can never guarantee the absence of high quality groups—which often is an insight equally
important as the presence of a strong pattern. For instance, in our example in Fig. 1,
it would be remarkable to establish that long-range interactions are to a large degree
independent of nanocluster geometry.
Therefore, in this paper (Sec. 3), we extend branch-and-bound search to objective
functions of the form
f (Q) = g(|Q|, med(Q), d(Q))
(2)
where g is monotonically increasing in the subpopulation size, monotonically decreasing
in any dispersion measure d around the median, and, besides that, depends only (but
in arbitrary form) on the subpopulation median. This involves developing an efficient
algorithm for computing the tight optimistic estimator given by the optimal value of the
objective function among all possible subsets of target values:
fˆ(Q) = max{f (R) : R ⊆ Q} ,

(3)

which has been shown to be a crucial ingredient for the practical applicability of branchand-bound (Grosskreutz et al, 2008; Lemmerich et al, 2016). So far, the most general
approach to this problem (first codified in Lemmerich et al (2016); generalized here in
Sec. 3.1) is to maintain a sorted list of target values throughout the search process and
then to compute Eq. (3) as the maximum of all subsets Ri ⊆ Q that contain all target
values of Q down to target value i—an algorithm that does not generalize to objective
functions depending on dispersion. This paper presents an alternative idea (Sec. 3.2)

3

where we do not fix the size of subset Ri as in the previous approach but instead fix
its median to target value i. It turns out that this suffices to efficiently compute the
tight optimistic estimator for all objective functions of the form of Eq. (2). Moreover,
we end up with a linear time algorithm (Sec. 3.3) in the important special case where
the dependence on size and dispersion is determined by the dispersion-corrected coverage
defined by


|Q|
amd(Q)
dcc(Q) =
max 1 −
,0
|P |
amd(P )
where amd denotes the mean absolute deviation from the median. This is the same computational complexity as the objective function itself. Consequently, this new approach
can discover subgroups according to a more refined selection criterion without increasing
the worst-case computational cost. Additionally, as demonstrated by empirical results
on a wide range of datasets (Sec. 4), it is also highly efficient and successfully reduces
the error of result subgroups in practice.

2. Subgroup Discovery
Before developing the novel approach to tight optimistic estimator computation, we
recall in this section the necessary basics of optimal subgroup discovery with numeric
target attributes. We focus on concepts that are essential from the optimization point
of view (see, e.g., Duivesteijn and Knobbe 2011 and references therein for statistical
considerations). As notional convention, we are using the symbol [m] for a positive
integer m to denote the set of integers {1, . . . , m}. Also, for a real-valued expression x
we write (x)+ to denote max{x, 0}. A summary of the most important notations used
in this paper can be found in Appendix C.

2.1. Description languages, objective functions, and closed selectors
Let P denote our given global population of entities, for each of which we know the
value of a real target variable y : P → R and additional descriptive information that
is captured in some abstract description language L of subgroup selectors σ : P →
{true, false}. Each of these selectors describes a subpopulation ext(σ) ⊆ P defined by
ext(σ) = {p ∈ P : σ(p) = true}
that is referred to as the extension of σ. Subgroup discovery is concerned with finding
descriptions σ ∈ L that have a useful (or interesting) distribution of target values in
their extension yσ = {y(p) : p ∈ ext(σ)}. This notion of usefulness is given by an
objective function f : L → R. That is, the formal goal is to find elements σ ∈ L
with maximal f (σ). Since we assume f to be a function of the multiset of y-values, let
us define f (σ) = f (ext(σ)) = f (yσ ) to be used interchangeably for convenience. One
example of a commonly used objective function is the impact measure ipa (see Webb
2001; here a scaled but order-equivalent version is given) defined by


mean(Q) − mean(P )
ipa(Q) = cov(Q)
(4)
max(P ) − mean(P ) +

4

where cov(Q) = |Q|/|P | denotes the coverage or relative size of Q (here—and wherever
else convenient—we identify a subpopulation Q ⊆ P with the multiset of its target
values).
The standard description language in the subgroup discovery literature1 is the language Lcnj consisting of logical conjunctions of a number of base propositions (or
predicates). That is, σ ∈ Lcnj are of the form
σ(·) ≡ πi1(·) ∧ · · · ∧ πil(·)
where the πij are taken from a pool of base propositions Π = {π1 , . . . , πk }. These
propositions usually correspond to equality or inequality constraints with respect to one
variable x out of a set of description variables {x1 , . . . , xn } that are observed for all
population members (e.g., π(p) ≡ x(p) ≥ v). However, for the scope of this paper it
is sufficient to simply regard them as abstract Boolean functions π : P → {true, false}.
In this paper, we focus in particular on the refined language of closed conjunctions
Ccnj ⊆ Lcnj (Pasquier et al, 1999), which is defined as Ccnj = {σ ∈ Lcnj : c(σ) = σ} by
the fixpoints of the closure operation c : Lcnj → Lcnj given by
^
c(σ) = {π ∈ Π : ext(π) ⊇ ext(σ)} .
(5)
These are selectors to which no further proposition can be added without reducing their
extension, and it can be shown that Ccnj contains at most one selector for each possible extension. While this can reduce the search space for finding optimal subgroups by
several orders of magnitude, closed conjunctions are the longest (and most redundant)
description for their extension and thus do not constitute intuitive descriptions by themselves. Hence, for reporting concrete selectors (as in Fig. 1), closed conjunctions have
to be simplified to selectors of approximately minimum length that describe the same
extension (Boley and Grosskreutz, 2009).

2.2. Branch-and-bound and optimistic estimators
The standard algorithmic approach for finding optimal subgroups with respect to a
given objective function is branch-and-bound search—a versatile algorithmic puzzle solving framework with several forms and flavors (see, e.g., Mehlhorn and Sanders, 2008,
Chap. 12.4). At its core, all of its variants assume the availability and efficient computability of two ingredients:
1. A refinement operator r : L → 2L that is monotone, i.e., for σ, ϕ ∈ L with
ϕ ∈ r(σ) it holds that ext(ϕ) ⊆ ext(σ), and that non-redundantly generates L.
That is, there is a root selector ⊥ ∈ L such that for every σ ∈ L there is a unique
sequence of selectors ⊥ = σ0 , σ1 , . . . , σl = σ with σi ∈ r(σi−1 ). In other words,
the refinement operator implicitly represents a directed tree (arborescence) on the
description language L rooted in ⊥.
1

In this article we remain with this basic setting for the sake of simplicity. It is, however, useful to
note that several generalizations of this concept have been proposed (e.g., Parthasarathy et al, 1999;
Huan et al, 2003), to which the contributions of this paper remain applicable.

5

Bst-BB(F, σ):
// F max. priority queue w.r.t. fˆ, σ current
f -maximizer
begin
if F = ∅ or fˆ(top(F))/f (σ) ≤ a then
return σ
else
R = r(top(F))
// refinement of fˆ-maximizer in queue
0
σ = argmax(f (ϕ) : ϕ ∈ {σ} ∪ R)
F 0 = (F \ {top(F)}) ∪ {ϕ ∈ R : fˆ(ϕ)/f (σ 0 ) ≥ a}
return Bst-BB(F 0 , σ 0 )
end
end
σ ∗ = Bst-BB({⊥}, ⊥)
// call with root element to find global
solution
Algorithm 1: Best-first branch-and-bound that finds a-approximation to objective function f based on refinement operator r and optimistic estimator fˆ; depthlimit and multiple solutions (top-k) parameters omitted; top denotes the find max
operation for priority queue.
2. An optimistic estimator (or bounding function) fˆ : L → R that bounds from
above the attainable subgroup value of a selector among all more specific selectors,
i.e., it holds that fˆ(σ) ≥ f (ϕ) for all ϕ ∈ L with ext(ϕ) ⊆ ext(σ).
Based on these ingredients, a branch-and-bound algorithm simply enumerates all elements of L starting from ⊥ using r (branch), but—based on fˆ—avoids expanding descriptions that cannot yield an improvement over the best subgroups found so far (bound).
Depending on the order in which language elements are expanded, one distinguishes between depth-first, breadth-first, breadth-first iterating deepening, and best-first search.
In the last variant, the optimistic estimator is not only used for pruning the search space,
but also to select the next element to be expanded, which is particularly appealing for
informed, i.e., tight, optimistic estimators. An important feature of branch-and-bound
is that it effortlessly allows to speed-up the search in a sound way by relaxing the result
requirement from being f -optimal to just being an a-approximation. That is, the
found solution σ satisfies for all σ 0 ∈ L that f (σ)/f (σ 0 ) ≥ a for some approximation
factor a ∈ (0, 1]. The pseudo-code given in Alg. 1 summarizes all of the above ideas.
Note that, for the sake of clarity, we omitted here some other common parameters such
as a depth-limit and multiple solutions (top-k), which are straightforward to incorporate
(see Lemmerich et al, 2016).
An efficiently computable refinement operator has to be constructed specifically for
the desired description language. For example for the language of conjunctions Lcnj , one
can define rcnj : Lcnj → Lcnj by
rcnj (σ) = {σ ∧ πi : max{j : πj ∈ σ} < i ≤ k}

6

where we identify a conjunction with the set of base propositions it contains. For the
closed conjunctions ccnj , let us define the lexicographical prefix of a conjunction σ ∈ Lcnj
and a base proposition index i ∈ [k] as σ|i = σ ∩ {π1 , . . . , πi }. Moreover, let us denote
with i(σ) the minimal index such that the i-prefix of σ is extension-preserving, i.e.,
i(σ) = min{i : ext(σ|i ) = ext(σ)}. With this we can construct a refinement operator
(Uno et al, 2004) rccj : Ccnj → 2Ccnj as
rccj (σ) = {ϕ : ϕ = ccnj (σ ∧ πj ), i(σ) < j ≤ k, πj 6∈ σ, ϕ|j = σ|j } .
That is, a selector ϕ is among the refinements of σ if ϕ can be generated by an application
of the closure operator given in Eq. (5) that is prefix-preserving.
How to obtain an optimistic estimator for an objective function of interest depends
on the definition of that objective. For instance, the coverage function cov is a valid
optimistic estimator for the impact function ipa as defined in Eq. (4), because the second
factor of the impact function is upper bounded by 1. In fact there are many different
optimistic estimators for a given objective function. Clearly, the smaller the value of
the bounding function for a candidate subpopulation, the higher is the potential for
pruning the corresponding branch from the enumeration tree. Ideally, one would like to
use fˆ(σ) = max{f (ϕ) : ext(ϕ) ⊆ ext(σ)}, which is the most strict function that still
is a valid optimistic estimator. In general, however, computing this function is as hard
as the original subgroup optimization problem we started with. Therefore, as a next
best option, one can disregard selectability and consider the (selection-unaware) tight
optimistic estimator given by
fˆ(σ) = max{f (R) : R ⊆ ext(σ)} .
This leaves us with a new combinatorial optimization problem: given a subpopulation
Q ⊆ P , find a sub-selection of Q that maximizes f . In the following section we will
discuss strategies for solving this optimization problem efficiently for different classes of
objective functions—including dispersion-corrected objectives.

3. Efficiently Computable Tight Optimistic Estimators
We are going to develop an efficient algorithm for the tight optimistic estimator in
three steps: First, we review and reformulate a general algorithm for the classic case
of non-dispersion-aware objective functions. Then we transfer the main idea of this
algorithm to the case of dispersion-corrected objectives based on the median, and finally
we consider a subclass of these functions where the approach can be computed in linear
time. Throughout this section we will identify a given subpopulation Q ⊆ P with the
multiset of its target values {y1 , . . . , ym } and assume that the target values are indexed
in ascending order, i.e., yi ≤ yj for i ≤ j. Also, it is helpful to define the following
partial order defined on finite multisets. Let Y = {y1 , . . . , ym } and Z = {z1 , . . . , zm0 } be
two multisets such that their elements are indexed in ascending order. We say that Y is
element-wise less or equal to Z and write Y ≤e Z if yi ≤ zi for all i ∈ [min{m, m0 }].

7

3.1. The standard case: monotone functions of a central tendency measure
The most general previous approach for computing the tight optimistic estimator for
subgroup discovery with a metric target variable is described by Lemmerich et al (2016)
where it is referred to as estimation by ordering. Here, we review this approach and give
a uniform and generalized version of that paper’s results. For this, we define the general
notion of a measure of central tendency as follows.
Definition 1. We call a mapping c : NR → R a (monotone) measure of central
tendency if for all multisets Y, Z ∈ NR with Y ≤e Z it holds that c(Y ) ≤ c(Z).
One can check that this definition applies to the standard measures of central tendency,
i.e., the arithmetic and geometric mean as well as the median2 med(Q) = ydm/2e , and
also to weighted variants of them (note, however, that it does not apply to the mode).
With this we can define the class of objective functions for which the tight optimistic
estimator can be computed efficiently by the standard approach as follows. We call
f : 2P → R a monotone level 1 objective function if it can be written as
f (Q) = g(|Q|, c(Q))
where c is some measure of central tendency and g is a function that is non-decreasing
in both of its arguments. One can check that the impact measure ipa falls under this
category of functions as do many of its variants.
The central observation for computing the tight optimistic estimator for monotone
level 1 functions is that the optimum value must be attained on a sub-multiset that
contains a consecutive segment of elements of Q from the top element w.r.t. y down to
some cut-off element. Formally, let us define the top sequence of sub-multisets of Q
as Ti = {ym−i+1 , . . . , ym } for i ∈ [m] and note the following observation:
Proposition 1. Let f be a monotone level 1 objective function. Then the tight optimistic
estimator of f can be computed as the maximum value on the top sequence, i.e., fˆ(Q) =
max{f (Ti ) : i ∈ [m]}.
Proof. Let R ⊆ Q be of size k with R = {yi1 , . . . , yik }. Since yij ≤ ym−j+1 , we have for
the top sequence element Tk that R ≤e Tk and, hence, c(R) ≤ c(Tk ) implying
f (R) = g(k, c(R)) ≤ g(k, c(Tl )) = f (Tk ) .
It follows that for each sub-multiset of Q there is a top sequence element of at least
equal objective value.
From this insight it is easy to derive an O(m) algorithm for computing the tight
optimistic estimator under the additional assumption that we can compute g and the
“incremental central tendency problem” (i, Q, (c(T1 ), . . . , c(Ti−1 )) 7→ c(Ti ) in constant
2

In this paper, we are using the simple definition of the median as the 0.5-quantile (as opposed to
defining it as (ym/2 + y1+m/2 )/2 for even m), which simplifies many of the definitions below and
additionally is well-defined in settings where averaging of target values is undesired.

8

time. Note that computing the incremental problem in constant time implies to only
access a constant number of target values and of the previously computed central tendency values. This can for instance be done for c = mean via the incremental formula
mean(Ti ) = ((i − 1) mean(Ti−1 ) + ym−i+1 )/i or for c = med through direct index access
of either of the two values ym−b(i−1)/2c or ym−d(i−1)/2e . Since, according to Prop. 1, we
have to evaluate f only for the m candidates Ti to find fˆ(Q) we can do so in time O(m)
by solving the problem incrementally for i = 1, . . . , m. The same overall approach can
be readily generalized for objective functions that are monotonically decreasing in the
central tendency or those that can be written as the maximum of one monotonically
increasing and one monotonically decreasing level 1 function. However, it breaks down
for objective functions that depend on more than just size and central tendency—which
inherently is the case when we want to incorporate dispersion-control.

3.2. Dispersion-corrected objective functions based on the median
We will now extend the previous recipe for computing the tight optimistic estimator to
objective functions that depend not only on subpopulation size and central tendency but
also on the target value dispersion in the subgroup. Specifically, we focus on the median as measure of central tendency and consider functions that are both monotonically
increasing in the described subpopulation size and monotonically decreasing in some dispersion measure around the median. To precisely describe this class of functions, we first
have to formalize the notion of dispersion measure around the median. For our purpose
the following definition suffices. Let us denote by Y∆med the multiset of absolute differences to the median of a multiset Y ∈ NR , i.e., Y∆med = {|y1 −med(Y )|, . . . , |ym −med(Y )|}.
Definition 2. We call a mapping d : NR → R a dispersion measure around the
median if d(Y ) is monotone with respect to the multiset of absolute differences to its
med then d(Y ) ≤ d(Z).
median Y∆med , i.e., if Y∆med ≤e Z∆
One can check that this definition contains the measures median absolute deviation
around the median mmd(Y ) = med(Y∆med ), the root mean of squared deviations around the
median rsm(Y ) = mean({x2 : x ∈ Y∆med })1/2 , as well as the mean absolute deviation
around the median amd(Y ) = mean(Y∆med ).3 Based on Def. 2 we can specify the class
of objective functions that we aim to tackle as follows: we call a function f : 2P → R a
dispersion-corrected or level 2 objective function (based on the median) if it can
be written as
f (Q) = g(|Q|, med(Q), d(Q))
(6)
where d is some dispersion measure around the median and g : R3 → R is a real function
that is non-decreasing in its first argument and non-increasing in its third argument
(without any monotinicity requirement for the second argument).
3

We work here with the given definition of dispersion measure because of its simplicity. Note, however,
that all subsequent arguments can be extended in a straightforward way to a wider class of dispersion
measures by considering the multisets of positive and negative deviations separately. This wider
class also contains the interquartile range and certain asymmetric measures, which are not covered
by Def. 2.

9

Q1

.
.
.
.
.
.
.

Q10

.
.
.
.
.
.
.
.
.

Q21

y1 · · · · · · · · · · · · · · · · · · · · y10 · · · · · · · · · · · · · · · · · · · · · · · · · · · · y17 · · · · · · · · · · · · · · · · · · · · · y21

Figure 2. Median sequence sets Q1 , . . . , Q21 (marked red) for 21 random target values y1 , . . . , y21
w.r.t. objective function f (Q) = |Q|/|P |−amd(Q)/amd(P ); the sets are identical for any arbitrary
dependence on the median med(Q) that could potentially be added to f , and for any such function
the optimal value is attained among those 21 sets (Prop. 2).

Our recipe for optimizing these functions is then to consider only subpopulations
R ⊆ Q that can be formed by selecting all individuals with a target value in some
interval. Formally, for a fixed index z ∈ {1, . . . , m} define mz ≤ m as the maximal
cardinality of a sub-multiset of the target values that has median index z, i.e.,
mz = min{2z, 2(m − z) + 1} .

(7)

Now, for k ∈ [mz ], let us define Qkz as the set with k consecutive elements around index
z. That is
n
o
Qkz = yz−b k−1 c , . . . , yz , . . . , yz+d k−1 e .
(8)
2

2

With this we can define the elements of the median sequence Qz as those subsets of∗
k
the form of Eq. (8) that maximize f for some fixed index z ∈ [m]. That is, Qz = Qz z
where kz∗ ∈ [mz ] is minimal with
∗

∗

f (Qkz z ) = g(kz∗ , yz , d(Qkz z )) = max{f (Qkz ) : k ∈ [mz ]} .
Thus, the number kz∗ is the smallest cardinality that maximizes the trade-off of size and
dispersion encoded by g (given the fixed median yz = med(Qkz ) for all k). Fig. 2 shows
an exemplary median sequence based on 21 random target values. In the following
proposition we note that, as desired, searching the median sequence is sufficient for
finding optimal subsets of Q.

10

Proposition 2. Let f be a dispersion-corrected objective function based on the median.
Then the tight optimistic estimator of f can be computed as the maximum value on the
median sequence, i.e., fˆ(Q) = max{f (Qz ) : z ∈ [m]}.
Proof. For a sub-multiset R ⊆ Q let us define the gap count γ(R) as
γ(R) = |{y ∈ Q \ R : min R < y < max R}| .
Let O ⊆ Q be an f -maximizer with minimal gap count, i.e., f (R) < f (O) for all R
with γ(R) < γ(O). Assume that γ(O) > 0. That means there is a y ∈ Q \ O such that
min O < y < max O. Define
(
(O \ {min O}) ∪ {y}, if y ≤ med(O)
S=
.
(O \ {max O}) ∪ {y}, otherwise
Per definition we have |S| = |O| and med(S) = med(O). Additionally, we can check that
med ≤ O med , and, hence, d(S) ≤ d(Q). This implies that
S∆
e ∆
f (S) = g(|S|, med(S), d(S)) ≥ g(|O|, med(O), d(O)) = f (O) .
However, per definition of S it also holds that γ(S) < γ(O), which contradicts that O
is an f -optimizer with minimal gap count. Hence, any f -maximizer O must have a gap
count of zero. In other words, O is of the form O = Qkz as in Eq. (8) for some median
z ∈ [m] and some cardinality k ∈ [mz ] and per definition we have f (Qz ) ≥ f (O) as
required.
Consequently, we can compute the tight optimistic estimator for any dispersioncorrected objective function based on the median in time O(m2 ) for subpopulations
of size m—again, given a suitable incremental formula for d. While this is not generally
a practical algorithm in itself, it is a useful departure point for designing one. In the
next section we show how it can be brought down to linear time when we introduce some
additional constraints on the objective function.

3.3. Reaching linear time—objectives based on dispersion-corrected
coverage
Equipped with the general concept of the median sequence, we can now address the
special case of dispersion-corrected objective functions where the trade-off between the
subpopulation size and target value dispersion is captured by a linear function of size
and the sum of absolute differences from the median. Concretely, let us define the
dispersion-corrected coverage (w.r.t. absolute median deviation) by




|Q|
amd(Q)
|Q| smd(Q)
dcc(Q) =
1−
=
−
|P |
amd(P ) +
|P | smd(P ) +
P
where smd(Q) = y∈Q |y − med(Q)| denotes the sum of absolute deviations from
the median. We then consider objective functions based on the dispersion-corrected
coverage of the form
f (Q) = g(dcc(Q), med(Q))
(9)

11

0.35

z =10
z =11
z =12
z =13

0.30

dcc( z

Qk )

0.25
0.20
0.15
0.10
0.05
0.00

5

k13∗ k11∗ k12∗
10

k

k10∗

15

20

Figure 3. Dispersion-corrected coverage of the sets Qkz as defined in Eq. (8) for median indices
z ∈ {10, 11, 12, 13} and the 21 random target values from Fig. 2; the sets Qz can be found in
∗
(Thm. 3).
incremental constant time since optimal size kz∗ is within a constant range of kz+1

where g is non-decreasing in its first argument. Let us note, however, that we could
replace the dcc function by any linear function that depends positively on |Q| and
negatively on smd. It is easy to verify that function of this form also obey the more
general definition of level-2 objective functions given in Sec. 3.2, and, hence can be
optimized via the median sequence.
The key to computing the tight optimistic estimator fˆ in linear time for functions
based on dispersion-corrected coverage is then that the members of the median sequence
Qz can be computed incrementally in constant time. Indeed, we can prove the following
theorem, which states that the optimal size for a multiset around median index z is
within 3 of the optimal size for a multiset around median index z + 1—a fact that can
also be observed in the example given in Fig. 2.
Theorem 3. Let f be of the form of Eq. (9). For z ∈ [m − 1] it holds for the size kz∗ of
the f -optimal multiset with median z that
∗
∗
kz∗ ∈ {max(0, kz+1
− 3), . . . , min(mz , kz+1
+ 3)} .

(10)

One idea to prove this theorem is to show that a) the gain in f for increasing the
multiset around a median index z is alternating between two discrete concave functions
and b) that the gains for growing multisets between two consecutive median indices are
bounding each other. For an intuitive understanding of this argument, Fig. 3 shows for
four different median indices z ∈ {10, 11, 12, 13} the dispersion-corrected coverage for the

12

let Q be given by {y1 , . . . , ym } in ascending order
compute el (i) and er (i) for i ∈ [m] through Eqs. (11) and (12)
∗ =1
f (Qm ) = g(1/|P |, ym ) and km
for z = m − 1 to 1 do
∗
∗
let k − = max(0, kz+1
− 3) and k + = min(mz , kz+1
+ 3) with mz as in Eq. (7)
−
+
for k = k to k do
let a = z − bk/2c and b = z + dk/2e
smd(Qkz ) = el (z) − el (a) − (a − 1)(yz − ya ) + er (z) − er (b) − (m − b)(yb − yz )
f (Qkz ) = g(k/|P | − smd(Qkz )/smd(P ), yz )
end
k∗

k∗

f (Qz ) = f (Qz z ) with kz∗ s.t. f (Qz z ) = max{f (Qkz ) : k − ≤ k ≤ k + }
end
fˆ(Q) = max{f (Qz ) : z ∈ [m]}
Algorithm 2: Linear time algorithm for computing tight optimistic estimator
fˆ(Q) of objective f (Q) = g(dcc(Q), med(Q)) as in Eq. (9).
sets Qkz as a function in k. On closer inspection, we can observe that when considering
only every second segment of each function graph, the corresponding dcc-values have
a concave shape. A detailed proof, which is rather long and partially technical, can be
found in Appendix A.
It follows that, after computing the objective value of Qm trivially as f (Qm ) =
g(1/|P |, ym ), we can obtain f (Qz−1 ) for z = m, . . . , 2 by checking the at most seven
candidate set sizes given by Eq. (10) as
o
n
kz+
kz−
z−1
f (Q ) = max f (Qz−1 ), . . . , f (Qz−1 )
with kz− = max(kz∗ − 3, 1) and kz+ = min(kz∗ + 3, mz ). It remains to see that we can compute individual evaluations of f in constant time (after some initial O(m) pre-processing
step). As a general data structure for quickly computing sums of absolute deviations
from a center point, we can define for i ∈ [m] the left error el (i) and the right error
er (i) as
i−1
m
X
X
el (i) =
yi − yj ,
er (i) =
yj − yi .
j=1

j=i+1

Note that we can compute these error terms for all i ∈ [m] in time O(m) via the
recursions
el (i) = el (i − 1) + (i − 1)(yi − yi−1 )

(11)

er (i) = er (i + 1) + (m − i)(yi+1 − yi )

(12)

and el (1) = er (m) = 0. Subsequently, we can compute sums of deviations from center
points of arbitrary subpopulations in constant time, as the following statement shows
(see Appendix B for a proof).

13

(med(Q)−med(P))/amd(P)

16
8

Q0
Q1

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

dataset id

Figure 4. Normalized median of optimal subgroup w.r.t. uncorrected positive median shift (Q0 )
and w.r.t. dispersion-corrected positive median shift (Q1 ) for 25 test datasets, sorted according
to median difference. Error bars show mean absolute median deviation of subgroups; groups
marked red have larger deviation than global deviation; fill color indicates group coverage from
0 (white) to 1 (black).

Proposition 4. Let Q = {y1 , . . . , ya , . . . , yz , . . . , yb , . . . , ym } be a multiset with 1 ≤ a <
z < b ≤ m and yi ≤ yj for i ≤ j. Then the sum of absolute deviations to yi of all
elements of the submultiset {ya , . . . , yz , . . . , yb } can be expressed as
b
X

|yz − yi | = el (z) − el (a) − (a − 1)(yz − ya )

i=a

+ er (z) − er (b) − (m − b)(yb − yz ) .
With this we can compute k 7→ f (Qkz ) in constant time (assuming g can be computed
in constant time). Together with Prop. 2 and Thm. 3 this results in a linear time
algorithm for computing Q 7→ fˆ(Q) (see Alg. 2 for a pseudo-code that summarizes all
ideas).

4. Dispersion-corrected Subgroup Discovery in Practice
The overall result of Sec. 3 is an efficient algorithm for dispersion-corrected subgroup
discovery which, e.g., allows us to replace the coverage term in standard objective functions by the dispersion-corrected coverage. To evaluate this efficiency claim as well
as the value of dispersion-correction, let us consider as objective the normalized and

14

15

abalone
ailerons
autoMPG8
baseball
california
compactiv
concrete
dee
delta ail
delta elv
elevators
forestfires
friedman
house
laser
mortgage
mv
pole
puma32h
stock
treasury
wankara
wizmir

Name

Target

delta E
Evdw-Evdw0

rings
goal
mpg
salary
med. h. value
usr
compr. strength
consume
sa
se
goal
area
output
price
output
30 y. rate
y
output
thetadd6
company10
1 m. def. rate
mean temp.
mean temp.
82 499
12, 200 250

69
357
24
24
72
202
70
60
66
66
155
70
48
160
42
128
79
260
318
80
128
87
82

|P | |Π|
4,177
13,750
392
337
20, 640
8,192
1,030
365
7, 129
9, 517
16, 599
517
1, 200
22, 784
993
1, 049
40, 768
14, 998
8, 192
950
1, 049
321
1, 461
0.106
0.131

9
−0.0008
22.5
740
179, 700
89
34.4
2.787
−0.0001
0.001
0.02
0.52
14.651
33, 200
46
6.71
-5.02086
0
0.000261
46.625
6.61
47.7
60

med(P )

0.277
0.088

2.359
0.000303
6.524
954.386
88, 354
9.661
13.427
0.831
0.000231
0.00198
0.00411
12.832
4.234
28,456
35.561
2.373
8.509
28.949
0.023
5.47
2.473
12.753
12.622
0.305
0.765

0.544
0.906
0.497
0.362
0.385
0.464
0.284
0.523
0.902
0.384
0.113
0.01
0.387
0.56
0.32
0.256
0.497
0.40
0.299
0.471
0.182
0.545
0.6

amd(Q1 )

aeff

|E0 |

0.378
0.34

0.43
0.217

0.202
0.234

0.373
0.081

0.118
0.0278

0.5 4, 712, 128
0.4 1, 498, 185

t0

t1

0.29
3.96

304
339
6, 542 460
0.11 0.09
0.22 0.21
2, 676 368
5, 161 208
43.9 35.8
2.05 1.29
33.3 6.11
8.9 4.01
13, 712 2, 891
23 23.7
3.23 1.56
7, 937 1, 308
0.96 0.83
11.6 1.59
31.9 13.2
2, 638
15
2, 648 15.5
12.5 1.56
14.8 5.91
11.9 1.24
38.5 1.48
204 1, 200
451 5, 650

690, 177
54, 103
67
117
65, 707
59, 053
221, 322
2, 653
2, 632
1, 415
526, 114
264, 207
2, 489
114, 566
815
1, 270
1, 017
2, 966
5, 782
1, 822
9, 247
4, 081
1, 409

|E1 |

Efficiency

0.191
11
11
2.257
1.662
1 848, 258
0.59 −0.0007 −0.0006 0.000288 0.000198 0.3 1, 069, 456
0.497
29
29
4.791
4.791
1
96
0.003
1550
2500 1245.092
0
1
117
0.019 262, 500 500, 001
94261 294, 00 0.4 1, 368, 662
0.603
94
93
7.8
3.472 0.5 2, 458, 105
0.1291
48.97
50.7 12.744
9.512
1 512, 195
0.381
3.815
4.008
0.721
0.434
1
18, 663
0.392 0.0001 0.0002 0.000226 0.000119
1
45, 194
0.369
0.002
0.002 0.00112 0.00108
1
10145
0.283
0.03
0.021 0.00813 0.00373 0.05 6, 356, 465
0.002
86.45 278.53 56.027
0
1 340, 426
0.294 18.934 19.727
3.065
2.73
1
19, 209
0.723 45, 200 34, 000 40, 576 27, 214 0.002 1, 221, 696
15.662
1
2, 008
0.093
109
135 40.313
0.097
11.61
14.41
2.081
0.98
1
40, 753
0.349
0.076
0.193
8.541
2.032
1
6, 513
0.24
100
100 38.995
16.692 0.2 1, 041, 146
0.244
0.026
0.031
0.018
0.017 0.4 3, 141, 046
0.337
52.5 54.375
3.741
2.515
1
85, 692
0.339 13.16
8.65
2.591
0.863
1
49, 197
0.296
60.6
67.6
8.873
4.752
1 191, 053
0.349
72.9
78.5
8.527
3.889
1 177, 768

amd(P ) cov(Q0 ) cov(Q1 ) med(Q0 ) med(Q1 ) amd(Q0 )

Selection Bias

Table 1. Datasets with corresponding population size (|P |), number of base propositions (|Π|), global median (med(P )) and mean absolute
median deviation (amd(P)) followed by coverage (cov(Q0 ), cov(Q1 )), median (med(Q0 ), med(Q1 )), and mean absolute median deviation
(amd(Q0 ), amd(Q1 )) for best subgroup w.r.t. non-dispersion corrected function f0 and dispersion-corrected function f1 , respectively; boldface indicates higher coverage and median and lower dispersion, underlines indicate higher dispersion than in global population; final
column segment contains accuracy parameter used in the efficiency study (aeff ) as well as number of expanded nodes (|E0 |, |E1 |) and
computation time in seconds (t0 , t1 ) for optimistic estimator based on top sequence fˆ0 and tight optimistic estimator fˆ1 , respectively—in
both cases when optimizing f1 ; depth-limit of 10 is used for all datasets with a < 1, no depth-limit otherwise.

24 binaries
25 gold

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

Dataset

dispersion-corrected impact function based on the median, i.e., f1 (Q) = dcc(Q)mds+ (Q)
where mds+ is the positive relative median shift


med(Q) − med(P )
mds+ (Q) =
.
max(P ) − med(P ) +
This function obeys Eq. (9); thus, its tight optimistic estimator can be computed using
the linear time algorithm from Sec. 3.3. The following empirical results were gathered by
applying it to a range of publicly available real-world datasets.4 We will first investigate
the effect of dispersion-correction on the output before turning to the effect of the tight
optimistic estimator on the computation time.

4.1. Selection Bias of Dispersion-Correction and its Statistical Merit
To investigate the selection bias of f1 let us also consider the non-dispersion corrected
variant f0 (Q) = cov(Q)mds+ (Q) where we simply replace the dispersion-corrected coverage by the ordinary coverage. This function is a monotone level 1 function, hence,
its tight optimistic estimator fˆ0 can be computed in linear time using the top sequence
approach. Fig. 4 shows the characteristics of the optimal subgroups that are discovered
with respect to both of these objective functions (see also Tab. 1 for exact values) where
for all datasets the language of closed conjunctions Ccnj has been used as description
language.
The first observation is that—as enforced by design—for all datasets the mean absolute
deviation from the median is lower for the dispersion-corrected variant (except in one
case where both functions yield the same subgroup). On average the dispersion for
f1 is 49 percent of the global dispersion, whereas it is 113 percent for f0 , i.e., when
not optimizing the dispersion it is on average higher in the subgroups than in the global
population. When it comes to the other subgroup characteristics, coverage and median
target value, the global picture is that f1 discovers somewhat more specific groups (mean
coverage 0.3 versus 0.44 for f0 ) with higher median shift (on average 0.73 normalized
median deviations higher). However, in contrast to dispersion, the behavior for median
shift and coverage varies across the datasets. In Fig. 4, the datasets are ordered according
to the difference in subgroup medians between the optimal subgroups w.r.t. f0 and those
w.r.t. f1 . This ordering reveals the following categorization of outcomes: When our
description language is not able to reduce the error of subgroups with very high median
value, f1 settles for more coherent groups with a less extreme but still outstanding central
tendency. On the other end of the scale, when no coherent groups with moderate size
and median shift can be identified, the dispersion-corrected objective selects very small
groups with the most extreme target values. The majority of datasets obey the global
trend of dispersion-correction leading to somewhat more specific subgroups with higher
median that are, as intended, more coherent.
4

Datasets contain all regression datasets from the KEEL repository (Alcalá et al, 2010) with at least
5 attributes and two materials datasets from the Nomad Repository nomad-coe.eu/; see Tab. 1.
Implementation available in open source Java library realKD bitbucket.org/realKD/. Computation
times determined on MacBook Pro 3.1 GHz Intel Core i7.

16

5

Pr(rope)

4

l̃1

3

24

2
1
4

0
1

1

11
14
12

0

1

2

l̃0

3

4

5

Pr(l0)

Pr(l1)

Figure 5. Effect of dispersion correction on lower bound of 95-percent confidence interval of
target variable; (left) improvement over global lower bound in standard deviations of dispersioncorrected objective (˜l1 ) and non-dispersion-corrected objective (˜l0 ) with annotations showing ids
of datasets where either method provides no improvement; (right) posterior joint probabilities
of the events that normalized difference (˜l1 − ˜l0 )/ max{˜l0 , ˜l1 } is larger than 0.1 (Pr(˜l1 )), less than
−0.1 (Pr(˜l0 )), or within [−0.1, 0.1] (Pr(rope)) according to Bayesian sign-test in barycentric
coordinates (sections correspond to regions where corresponding event is maximum a posteriori
outcome).

To determine based on these empirical observations, whether we should generally favor
dispersion correction, we have to specify an application context that specifies the relative
importance of coverage, central tendency, and dispersion. For that let us consider the
common statistical setting in which we do not observe the full global population P
but instead subgroup discovery is performed only on an i.i.d. sample P 0 ⊆ P yielding
subpopulations Q0 = σ(P 0 ). While σ has been optimized w.r.t. the statistics on that
sample Q0 we are actually interested in the properties of the full subpopulation Q = σ(P ).
For instance, a natural question is what is the minimal y-value that we expect to see
in a random individual q ∈ Q with high confidence. That is, we prefer subgroups with
an as high as possible threshold l such that a random q ∈ Q satisfies with probability5
1 − δ that y(q) ≥ l. This criterion gives rise to a natural trade-off between the three
evaluation metrics through the empirical Chebycheff inequality (see Kabán, 2012,
Eq. (17)), according to which we can compute such a value as mean(Q0 ) − (Q0 ) where
s
(|Q0 |2 − 1)var(Q0 )
(Q0 ) =
|Q0 |2 δ − |Q0 |
P
and var(Y ) = y∈Y (y − mean(Y ))2 /(|Y | − 1) is the sample variance. Note that this
expression is only defined for sample subpopulations with a size of at least 1/δ. For
5

The probability is w.r.t. to the distribution with which the sample P 0 ⊆ P is drawn.

17

smaller subgroups our best guess for a threshold value would be the one derived from
the global sample mean(P 0 ) − (P 0 ) (which we assume to be large enough to determine
an -value). This gives rise to the following standardized lower confidence bound
score ˜l that evaluates how much a subgroup improves over the global l value:
!
(
0 ) − l(P 0 )
mean(Q0 ) − (Q0 ) , if (Q0 ) defined
l(Q
˜l(Q0 ) =
p
where l(Q0 ) =
.
mean(P 0 ) − (P 0 ) , otherwise
var(P 0 )
+

The plot on the left side of Fig. 5 shows the score values of the optimal subgroup
w.r.t. to f1 (˜l1 ) and f0 (˜l0 ) using confidence parameter δ = 0.05. Except for three
exceptions (datasets 3,4, and 12), the subgroup resulting from f1 provides a higher
lower bound than those from the non-dispersion corrected variant f0 . That is, the data
shows a strong advantage for dispersion correction when we are interested in selectors
that mostly select individuals with a high target value from the underlying population P .
In order to test the significance of these results, we can employ the Bayesian sign-test
(Benavoli et al, 2014), a modern alternative to classic frequentist null hypothesis tests
that avoids many of the well-known disadvantages of those (see Demšar, 2008; Benavoli
et al, 2016). With Bayesian hypothesis tests, we can directly evaluate the posterior
probabilities of hypotheses given our experimental data instead of just rejecting a null
hypothesis based on some arbitrary significance level. Moreover, we differentiate between
sample size and effect size by the introduction of a region of practical equivalence (rope).
Here, we are interested in the relative difference z̃ = (˜l1 − ˜l0 )/(max{˜l0 , ˜l1 }) on average for
random subgroup discovery problems. Using a conservative choice for the rope, we call
the two objective functions practically equivalent if the mean z̃-value is at most r = 0.1.
Choosing the prior belief that f0 is superior, i.e., z̃ < −r, with a prior weight of 1, the
procedure yields based on our 25 test datasets the posterior probability of approximately
1 that z̃ > r on average (see the right part of Fig. 5 for in illustration of the posterior
belief). Hence, we can conclude that dispersion-correction improves the relative lower
confidence bound of target values on average by more than 10 percent when compared
to the non-dispersion-corrected function.

4.2. Efficiency of the Tight Optimistic Estimator
To study the effect of the tight optimistic estimator, let us compare its performance
to that of a baseline estimator that can be computed with the standard top sequence
approach. Since f1 is upper bounded by f0 , fˆ0 is a valid, albeit non-tight, optimistic
estimator for f1 and can thus be used for this purpose. The exact speed-up factor
is determined by the ratio of enumerated nodes for both variants as well as the ratio
of computation times for an individual optimistic estimator computation. While both
factors determine the practically relevant outcome, the number of nodes evaluated is
a much more stable quantity, which indicates the full underlying speed-up potential
independent of implementation details. Similarly, “number of nodes evaluated” is also
an insightful unit of time for measuring optimization progress. Therefore, in addition to
the computation time in seconds t0 and t1 , let us denote by E0 , E1 ⊆ L the set of nodes

18

105

fˆ0 (σ0+(n))
fˆ1 (σ1+(n))

t0 /t1
|E0 |/|E1 |

f1 (σ1∗(n))
f1 (σ0∗(n))

103

value

speed-up

104

10-1

102
101
100

10-2 4
10

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

10-1

dataset id

105

n (#nodes)

106

Figure 6. Effect of tight optimistic estimator; (right) optimistic estimation (fˆ1 (σ1+ ), fˆ0 (σ0+ ))
of remaining search space and value of current best solution (f1 (σ1∗ ), f1 (σ0∗ )) resulting from tight
estimator and top sequence estimator, respectively—per processed nodes for dataset binaries;
(left) speedup factor (t0 /t1 ) in order for all datasets plus potential speed-up factor (|E0 |/|E1 |).

enumerated by branch-and-bound using fˆ0 and fˆ1 , respectively—but in both cases for
optimizing the dispersion-corrected objective f1 . Moreover, when running branch-andbound with optimistic estimator fˆi , let us denote by σi∗ (n) and σi+ (n) the best selector
found and the top element of the priority queue (w.r.t. fˆi ), respectively, after n nodes
have been enumerated.
The plot on the left side of Fig. 6 shows the speed-up factor t1 /t0 on a logarithmic
axis for all datasets in increasing order along with the potential speed-up factors |E0 |/|E1 |
(see Tab. 1 for numerical values). There are seven datasets for which the speed-up turns
out to be minor followed by four datasets with a modest speed-up factor of 2. For the
remaining 14 datasets, however, we have solid speed-up factors between 4 and 20 and
in four cases immense values between 100 and 4, 000. This demonstrates the decisive
potential effect of tight value estimation even when compared to another non-trivial
estimator like fˆ0 (which itself improves over simpler options by orders of magnitude;
see Lemmerich et al 2016). Similar to the results in Sec. 4.1, the Bayesian sign-test
for the normalized difference z = (t1 − t0 )/ max{t1 , t0 } with the prior set to practical
equivalence (z ∈ [−0.1, 0.1]) reveals that the posterior probability of fˆ1 being superior
to fˆ0 is approximately 1.
In almost all cases the potential speed-up given by the ratio of enumerated nodes
is considerably higher than the actual speed-up, which shows that, despite the same
asymptotic time complexity, an individual computation of the tight optimistic estimator
is slower than the simpler top sequence based estimator—but also indicates that there
is room for improvements in the implementation. When zooming in on the optimization
progress over time for the binaries dataset, which exhibits the most extreme speed-up

19

(right plot in Fig. 6), we can see that not only does the tight optimistic estimator close the
gap between best current selector and current highest potential selector much faster—
thus creating the huge speed-up factor—but also that it causes better solutions to be
found earlier. This is an important property when we want to use the algorithm as an
anytime algorithm, i.e., when allowing the user to terminate computation preemptively,
which is important in interactive data analysis systems. This is an advantage enabled
specifically by using the tight optimistic estimators in conjunction with the best-first
node expansion strategy.

5. Conclusion
During the preceding sections, we developed and evaluated an effective algorithm for
simultaneously optimizing size, central tendency, and dispersion in subgroup discovery
with a numerical target. This algorithm is based on two central results: 1) the tight
optimistic estimator for any objective function that is based on some dispersion measure around the median can be computed as the function’s maximum on a linear-sized
sequence of sets—the median sequence (Prop. 2); and 2) for objective functions based
on the concept of the dispersion-corrected coverage w.r.t. the absolute deviation from
the median, the individual sets of the median sequence can be generated in incremental
constant time (Thm. 3).
Among the possible applications of the proposed approach, the perhaps most important one is to replace the standard coverage term in classic objective functions by the
dispersion-corrected coverage, i.e., the relative subgroup size minus the relative subgroup
dispersion, to reduce the error of result subgroups—where error refers to the descriptive
or predictive inaccuracy incurred when assuming the median value of a subgroup for
all its members. As we saw empirically for the impact function (based on the median),
this correction also has a statistical advantage resulting in subgroups where we can assume greater target values for unseen group members with high confidence. In addition
to enabling dispersion-correction to known objective functions, the presented algorithm
also provides novel degrees of freedom, which might be interesting to exploit in their
own right: The dependence on the median is not required to be monotone, which allows
to incorporate a more sophisticated influence of the central tendency value than simple
monotone average shifts. For instance, given a suitable statistical model for the global
distribution, the effect of the median could be a function of the probability P[med(Q)],
e.g., its Shannon information content. Furthermore, the feasible dispersion measures
allow for interesting weighting schemes, which include possibilities of asymmetric effects
of the error (e.g., for only punishing one-sided deviation from the median).
Regarding the limitations of the presented approach, let us note that it cannot be directly applied
p to the previously proposed dispersion-aware functions, i.e., the t-score
tsc(Q) =
|Q|(mean(Q) − mean(P ))/std(Q) and the mmad score for ranked data
mmd(Q) = |Q|/(2med(Q) + mmd(Q)). While both of these functions can be optimized
via the median sequence approach (assuming a t-score variant based on the median),
we are lacking an efficient incremental formula for computing the individual function

20

values for all median sequence sets, i.e., a replacement for Thm. 3. Though finding such
a replacement in future research is conceivable, this leaves us for the moment with a
quadratic time algorithm (in the subgroup size) for the tight optimistic estimator, which
is not generally feasible (although potentially useful for smaller datasets or as part of
a hybrid optimistic estimator, which uses the approach for sufficiently small subgroups
only).
Since they share basic monotonicities, it is possible to use functions based on dispersioncorrected coverage as an optimization proxy for the above mentioned objectives. For
instance, the ranking of the top
p 20 subgroups w.r.t. the dispersion-corrected binomial
quality function, dcb(Q) =
dcc(Q)(med(Q) − med(P )), turns out to have a mean
Spearman rank correlation coefficient with the median-based t-score of apx. 0.783 on
five randomly selected test datasets (delta elv, laser, stock, treasury, gold ). However, a
more systematic understanding of the differences and commonalities of these functions is
necessary to reliably replace them with one another. Moreover, the correlation deteriorates quite sharply when we compare to the original mean/variance based t-score (mean
Spearman correlation coefficient 0.567), which points to the perhaps more fundamental
limitation of the presented approach for dispersion-correction: it relies on using the median as measure of central tendency. While the median and the mean absolute deviation
from the median are an interpretable, robust, and sound combination of measures (the
median of a set of values minimizes the sum of absolute deviations), the mean and the
variance are just as sound, are potentially more relevant when sensitivity to outliers
is required, and provide a wealth of statistical tools (e.g., the empirical Chebyshev’s
inequality used above).
Hence, a straightforward but valuable direction for future work is the extension of efficient tight optimistic estimator computation to dispersion-correction based on the mean
and variance. A basic observation for this task is that objective functions based on dispersion measures around the mean must also attain their maximum on gap-free intervals
of target values. However, for a given collection of target values, there is a quadratic
number of intervals such that a further idea is required in order to attain an efficient,
i.e., (log-)linear time algorithm. Another valuable direction for future research is the
extension of consistency and error optimization to the case of multidimensional target
variables where subgroup parameters can represent complex statistical models (known
as exceptional model mining Duivesteijn et al, 2016). While this setting is algorithmically more challenging than the univariate case covered here, the underlying motivation
remains: balancing group size and exceptionality, i.e., distance of local to global model
parameters, with consistency, i.e., local model fit, should lead to the discovery of more
meaningful statements about the data and the underlying domain.
Acknowledgements The authors thank the anonymous reviewers for their useful and
constructive suggestions. Jilles Vreeken and Mario Boley are supported by the Cluster
of Excellence “Multimodal Computing and Interaction” within the Excellence Initiative
of the German Federal Government. Bryan R. Goldsmith acknowledges support from
the Alexander von Humboldt-Foundation with a Postdoctoral Fellowship. Additionally,

21

this work was supported through the European Union’s Horizon 2020 research and innovation program under grant agreement no. 676580 with The Novel Materials Discovery
(NOMAD) Laboratory, a European Center of Excellence.

References
Alcalá J, Fernández A, Luengo J, Derrac J, Garcı́a S, Sánchez L, Herrera F (2010) Keel
data-mining software tool: Data set repository, integration of algorithms and experimental analysis framework. Journal of Multiple-Valued Logic and Soft Computing
17(2-3):255–287
Atzmueller M (2015) Subgroup discovery. Wiley Interdisciplinary Reviews: Data Mining
and Knowledge Discovery 5(1):35–49
Bay SD, Pazzani MJ (2001) Detecting group differences: Mining contrast sets. Data
Mining and Knowledge Discovery 5(3):213–246
Benavoli A, Corani G, Mangili F, Zaffalon M, Ruggeri F (2014) A bayesian wilcoxon
signed-rank test based on the dirichlet process. In: ICML, pp 1026–1034
Benavoli A, Corani G, Demsar J, Zaffalon M (2016) Time for a change: a tutorial for comparing multiple classifiers through bayesian analysis. arXiv preprint arXiv:160604316
Boley M, Grosskreutz H (2009) Non-redundant subgroup discovery using a closure system. In: Joint European Conf. on Machine Learning and Knowledge Discovery in
Databases, Springer, pp 179–194
Boley M, Moens S, Gärtner T (2012) Linear space direct pattern sampling using coupling
from the past. In: Proc. of the 18th ACM SIGKDD int. conf. on Knowledge discovery
and data mining, ACM, pp 69–77
Demšar J (2008) On the appropriateness of statistical tests in machine learning. In:
Workshop on Evaluation Methods for Machine Learning in conjunction with ICML
Duivesteijn W, Knobbe A (2011) Exploiting false discoveries–statistical validation of
patterns and quality measures in subgroup discovery. In: IEEE 11th Int. Conf. on
Data Mining, IEEE, pp 151–160
Duivesteijn W, Feelders AJ, Knobbe A (2016) Exceptional model mining. Data Mining
and Knowledge Discovery 30(1):47–98
Friedman JH, Fisher NI (1999) Bump hunting in high-dimensional data. Statistics and
Computing 9(2):123–143
Goldsmith BR, Boley M, Vreeken J, Scheffler M, Ghiringhelli LM (2017) Uncovering
structure-property relationships of materials by subgroup discovery. New Journal of
Physics 19(1):13–31

22

Grosskreutz H, Rüping S, Wrobel S (2008) Tight optimistic estimates for fast subgroup
discovery. In: Joint European Conf. on Machine Learning and Knowledge Discovery
in Databases, Springer, pp 440–456
Grosskreutz H, Boley M, Krause-Traudes M (2010) Subgroup discovery for election analysis: a case study in descriptive data mining. In: Int. Conf. on Discovery Science,
Springer, pp 57–71
Huan J, Wang W, Prins J (2003) Efficient mining of frequent subgraphs in the presence
of isomorphism. In: 3rd IEEE Int. Conf. on Data Mining, IEEE, pp 549–552
Kabán A (2012) Non-parametric detection of meaningless distances in high dimensional
data. Statistics and Computing 22(2):375–385
Klösgen W (1996) Explora: A multipattern and multistrategy discovery assistant. In:
Advances in knowledge discovery and data mining, American Association for Artificial
Intelligence, pp 249–271
Klösgen W (2002) Data mining tasks and methods: Subgroup discovery: deviation
analysis. In: Handbook of data mining and knowledge discovery, Oxford University
Press, Inc., pp 354–361
Lavrač N, Kavšek B, Flach P, Todorovski L (2004) Subgroup discovery with cn2-sd.
Journal of Machine Learning Research 5(Feb):153–188
Lemmerich F, Atzmueller M, Puppe F (2016) Fast exhaustive subgroup discovery with
numerical target concepts. Data Mining and Knowledge Discovery 30(3):711–762
Li G, Zaki MJ (2016) Sampling frequent and minimal boolean patterns: theory and
application in classification. Data Mining and Knowledge Discovery 30(1):181–225
Mehlhorn K, Sanders P (2008) Algorithms and data structures: The basic toolbox.
Springer Science & Business Media
Parthasarathy S, Zaki MJ, Ogihara M, Dwarkadas S (1999) Incremental and interactive
sequence mining. In: Proc. 8th int. conf. on Information and knowledge management,
ACM, pp 251–258
Pasquier N, Bastide Y, Taouil R, Lakhal L (1999) Efficient mining of association rules
using closed itemset lattices. Information systems 24(1):25–46
Pieters BF, Knobbe A, Dzeroski S (2010) Subgroup discovery in ranked data, with an
application to gene set enrichment. In: Proc. preference learning workshop (PL 2010)
at ECML PKDD, vol 10, pp 1–18
Schmidt J, Hapfelmeier A, Mueller M, Perneczky R, Kurz A, Drzezga A, Kramer S
(2010) Interpreting pet scans by structured patient data: a data mining case study in
dementia research. Knowledge and Information Systems 24(1):149–170

23

Song H, Kull M, Flach P, Kalogridis G (2016) Subgroup discovery with proper scoring
rules. In: Joint European Conference on Machine Learning and Knowledge Discovery
in Databases, Springer, pp 492–510
Uno T, Asai T, Uchida Y, Arimura H (2004) An efficient algorithm for enumerating
closed patterns in transaction databases. In: Int. Conf. on Discovery Science, Springer,
pp 16–31
Webb GI (1995) Opus: An efficient admissible algorithm for unordered search. Journal
of Artificial Intelligence Research 3:431–465
Webb GI (2001) Discovering associations with numeric variables. In: Proc. of the 7th
ACM SIGKDD int. conf. on Knowledge discovery and data mining, ACM, pp 383–388
Wrobel S (1997) An algorithm for multi-relational discovery of subgroups. In: European
Symposium on Principles of Data Mining and Knowledge Discovery, Springer, pp
78–87

A. Proof of Theorem 3
In order to proof Thm. 3, let us start by noting that for functions of the form of
Eq. (9), finding the set size kz∗ corresponds to maximizing the dispersion-corrected coverage among all multisets with consecutive elements around median yz (as defined in
Eq. 8). In order to analyze this problem, let us write
hz (k) = dcc(Qkz ) =

|Qkz | smd(Qkz )
−
|P |
smd(P )

for the dispersion-corrected coverage of the multiset Qkz . Let ∆hz : [mz ] → R denote
the difference or gain function of hz , i.e., ∆hz (k) = hz (k) − hz (k − 1) where we consider
Q0z = ∅ and, hence, hz (0) = 0. With this definition we can show that hz is alternating
between two concave functions, i.e., considering either only the even or only the odd
subset of its domain, the gains are monotonically decreasing. More precisely:
Lemma 5. For all k ∈ [mz ] \ {1, 2} we have that ∆hz (k) ≤ ∆hz (k − 2).
Proof. For k ∈ [mz ], let us denote by qzk the additional y-value that Qkz contains compared
to Qk−1
(considering Q0z = ∅), i.e., Qkz \ Qk−1
= {qzk }. We can check that
z
z

q
, k odd
z−b k−1
2 c
qzk =
.
qz+ k−1 , k even
d 2 e

24

With this and using the shorthands n = |P | and d = smd(P ) we can write
∆hz (k) − ∆hz (k − 2) = hz (k) − hz (k − 1) − (hz (k − 2) − hz (k − 3))
k smd(Qkz ) k − 1 smd(Qk−1
) k − 2 smd(Qk−2
) k − 3 smd(Qk−3
)
z
z
z
−
−
+
−
+
+
−
n
d
n
d
n
d
n
d

1
1
= (k − k + 1 − k + 2 + k − 3) +
smd(Qk−2
) − smd(Qkz ) + smd(Qk−1
) − smd(Qk−3
)
z
z
z
{z
} d
n|
=

=0


1
=
−|qzk − yz | − |qzk−1 − yz | + |qzk−1 − yz | + |qzk−2 − yz |
d

1
=
−|qzk − yz | + |qzk−2 − yz |
d
case k odd
 

1 
=
− yz − yz−b k−1 c + yz − yz−b k−3 c
= yz−b k−1 c − yz−b k−3 c ≤ 0
2
2
2
2
d
case k even
 

1 
=
− yz+d k−1 e − yz + yz+d k−3 e − yz
= yz+d k−3 e − yz+d k−1 e ≤ 0
2
2
2
2
d

One important consequence of this fact is that the operation of growing a set around
median z by two elements—one to the left and one to the right—has monotonically
decreasing gains. In other words, the smoothed function hz (k) = hz (k) + hz (k − 1) is
concave or formally
∆hz (k) + ∆hz (k − 1) ≥ ∆hz (k + 1) + ∆hz (k) .

(13)

Moreover, we can relate the gain functions of consecutive median indices as follows.
Lemma 6. Let z ∈ [m] \ {1} and k ∈ [mz−1 ] \ {1, 2, 3}. It holds that
∆hz−1 (k − 2) + ∆hz−1 (k − 3) ≥ ∆hz (k) + ∆hz (k − 1)
∆hz−1 (k) + ∆hz−1 (k − 1) ≤ ∆hz (k − 2) + ∆hz (k − 3)

(14)
(15)

Proof. For this proof, let us use the same shorthands as in the proof of Lemma 5 and
start by noting that for all i ∈ [m] and k ∈ [mz ] \ {1} we have the equality
∆hi (k) + ∆hi (k − 1) =

2 |qik − yi | + |qik−1 − yi |
−
n
d

(16)

which we can see by extending
∆hi (k) + ∆hi (k − 1) = hi (k) − hi (k − 1) + hi (k − 1) − hi (k − 2)
=

k − k + 2 smd(Qki ) − smd(Qik−2 )
2 |q k − yi | + |qik−1 − yi |
−
= − i
.
n
d
n
d

25

We can then show Eq. (14) by applying Eq. (16) two times to
∆hz−1 (k − 2) + ∆hz−1 (k − 3) − (∆hz (k) + ∆hz (k − 1))

1
k−2
k−3
−|qz−1
− yz−1 | − |qz−1
− yz−1 | + |qzk − yz | + |qzk−1 − yz |
=
d
and finally by checking separately the case k odd
=


1
yz−1−b k−3 c − yz−1 + yz−1 − yz−1+d k−4 e + yz − yz−b k−1 c + yz+d k−2 e − yz
2
2
2
2
d



1

= 
−
y
y
k−1 + y
k−1 − y
k
k
z−1+
z−3+
z−
z−
d2e
d 2 e ≥ 0
b 2 c
d | b 2 c {z
{z
}
} |
≥0

=0

and the case k even
=


1
yz−1 − yz−1+d k−3 e + yz−1−b k−4 c − yz−1 + yz+d k−1 e − yz + yz − yz−b k−2 c
2
2
2
2
d



1

y
= 
−
y
k−1 − y
k−1
k
k +y
z+1−
z+1−
z+
z−2+
b2c
b2c
d 2 e
d 2 e ≥ 0 .
d |
{z
} |
{z
}
≥0

=0

Similarly, for Eq. (15) by applying Eq. (16) two times we can write
∆hz−1 (k) + ∆hz−1 (k − 1) − (∆hz (k − 2) + ∆hz (k − 3))

1
k−1
k
− yz−1 | + |qzk−2 − yz | + |qzk−3 − yz |
=
−|qz−1
− yz−1 | − |qz−1
d 








1


y
−
y
k−1 − y
k−1 + y
k
k

z−2+
z−2+
z−1−
z+1−
d


d 2 e  ≤ 0, k odd
d2e
b 2 c
b 2 c


|
{z
}
|
{z
}

=0
≤0


=







1


y
−
y
−
y
k even
k−1
k−1
k
k +y

z+2−b 2 c
z−1+d 2 e  ≤ 0,
z−1+d 2 e
 d  z−b 2 c


|
{z
} |
{z
}

≤0

=0

Combining all of the above we can finally proof our main result as follows.
Theorem 3. We start by showing that every k ∈ [mz+1 ] with k < kz∗ − 3 can not be an
∗ , and, hence, k ∗ ≤ k ∗
optimizer of hz+1 . It follows that kz∗ − 3 ≤ kz+1
z
z+1 + 3 as required

26

for the upper bound. Indeed, we have
hz+1 (k) = hz+1 (k + 2) − (∆hz+1 (k + 2) + ∆hz+1 (k + 1))
≤ hz+1 (k + 2) − (∆hz+1 (kz∗ − 2) + ∆hz+1 (kz∗ − 3))
≤ hz+1 (k + 2) −

(∆hz (kz∗ )
|

+

∆hz (kz∗
{z

>0 by def. of

kz∗

(by Eq. (13))

− 1)) < hz+1 (k + 2) .
}

(by Lm. 6)

Analogously, for the lower bound, we show that every k ∈ [mz+1 ] with k > kz∗ + 3
∗ , and, hence,
can not be the smallest optimizer of hz+1 . It follows that kz∗ + 3 ≥ kz+1
∗
kz∗ ≥ kz+1
− 3 as required. Indeed, we can write
hz+1 (k) = hz+1 (k − 2) + ∆hz+1 (k) + ∆hk+1 (k − 1)
≤ hz+1 (k − 2) + ∆hz+1 (kz∗ + 4) + ∆hz+1 (kz∗ + 3)
≤ hz+1 (k − 2) +

∆hz (kz∗
|

+ 2) +
{z

∆hz (kz∗

≤0 by def. of kz∗

(by Eq. (13))

+ 1) ≤ hz+1 (k − 2)
}

(by Lm. 6)

B. Additional Proofs
Prop. 4. Using dij as a shorthand for yj − yi for i, j ∈ [m] with i ≤ j we can write
el (z) − el (a) − (a − 1)(yz − ya ) + er (z) − er (b) − (m − b)dzb
=

=

z−1
X

diz −

i=1

i=1

z−1
X

a−1
X

diz +

i=a

=

a−1
X

z−1
X

i=1

m
X

dia − (a − 1)daz +

(d − d ) −(a − 1)daz +
| iz {z ia}

diz + (a − 1)daz − (a − 1)daz +

i=a

=

z−1
X
i=a

dzi −

i=z+1

daz

b
X
i=z+1

dzi =

b
X

dbi − (m − b)dzb

i=b+1
b
X
i=z+1

b
X
i=z+1

diz +

m
X

|yz − yi |

i=a

27

dzi +

m
X
i=b+1

(dzi − dbi ) −(m − b)dzb
| {z }
dzb

dzi + (m − b)dzb − (m − b)dzb

C. Summary of Used Notations
Symbol

Meaning

Defined in

|·|
[k]
(x)+
≤e
2X
NX
σ, ϕ
c
d
el (i), er (i)
f
fˆ

cardinality of a set or absolute value of a number
set of integers {1, . . . , k}
max{x, 0} for a real-valued expression x
element-wise less-or-equal relation for multisets of real values
power set of a set X, i.e., set of all of its subsets
set of all multisets containing elements from set X
subgroup selectors σ, ϕ : P → {true, false}
measure of central tendency
measure of dispersion
left and right cumulative errors of target values up to value i
objective function
tight optimistic estimator of objective function f
number of elements in subpopulation Q
maximal size parameter k for consecutive value set Qkz
f -maximizing size parameter k for consecutive value set Qkz
numeric target attribute y : P → R
i-th target value of subpopulation w.r.t. ascending order
global population of given subgroup discovery problem
some subpopulation Q ⊆ P
median sequence element with median index z
submultiset of Q with k consecutive elements around index z
top sequence element i, i.e., Ti = {ym−i+1 , . . . , ym }
real-valued multiset
multiset of differences of elements in Y to its median
description language and language of conjunctions
language of closed conjunctions
mean absolute deviation of y-values in Q to their median
coverage, i.e., relative size |Q|/|P | of subpopulation Q
dispersion-corrected coverage of subpopulation Q
arithmetic mean of y-values in Q
impact, i.e., weighted mean-shift, of subpopulation Q
median of y-values in Q
sum of absolute deviations of y-values in Q to their median

2
2
3
2.1
3.1
3.2
3.3
2.1
2.2
3
3.2
3.2
2.1
3
2.1
2.1
3.2
3.2
3.1
3
3.2
2.1
2.1
3.2
2.1
3.3
2.1
3.1
3.3

m
mz
kz∗
y
yi
P
Q
Qz
Qkz
Ti
Y
Y∆med
L, Lcnj
Ccnj
amd(Q)
cov(Q)
dcc(Q)
mean(Q)
ipa(Q)
med(Q)
smd(Q)

28

