Quantifying Program Bias
Aws Albarghouthi, Loris D’Antoni, and
Samuel Drews

Aditya Nori
Microsoft Research

arXiv:1702.05437v2 [cs.PL] 7 Mar 2017

University of Wisconsin–Madison

Abstract
With the range and sensitivity of algorithmic decisions expanding at a break-neck speed, it is imperative that we
aggressively investigate whether programs are biased. We propose a novel probabilistic program analysis technique and
apply it to quantifying bias in decision-making programs. Specifically, we (i) present a sound and complete automated
verification technique for proving quantitative properties of probabilistic programs; (ii) show that certain notions of
bias, recently proposed in the fairness literature, can be phrased as quantitative correctness properties; and (iii) present
FairSquare, the first verification tool for quantifying program bias, and evaluate it on a range of decision-making programs.

1.

Introduction

A number of very interesting applications of program analysis have been explored in the probabilistic setting: reasoning
about cyber-physical systems [66], proving differential privacy of complex algorithms [14], reasoning about approximate
programs and hardware [23, 65], synthesizing control programs [25], amongst many others. In this paper, we turn our
attention to the problem of quantifying program bias.
Program bias Programs that make decisions can be biased. Consider, for instance, automatic grading of writing prompts
for standardized tests [10]; some speech patterns may be characterized as poor writing style and result in lower scores.
However, if such speech patterns are affiliated with a specific ethnic group, then the bias is a potential source of concern.
Programs have become powerful arbitrators of a range of significant decisions with far-reaching societal impact—
hiring [50, 55], welfare allocation [40], prison sentencing [8], policing [19, 62], amongst many others. With the range and
sensitivity of algorithmic decisions expanding by the day, the problem of understanding the nature of program bias is a
pressing one: Indeed, the notion of algorithmic fairness has recently captured the attention of a broad spectrum of experts,
within computer science and without [7, 8, 12, 22, 32, 37, 41, 67–69, 72].
Fairness and justice have always been a ripe topic for philosophical debate [63], and, of course, there are no established
rigorous definitions. Nonetheless, the rise of automated decision-making prompted the introduction of a number of formal
definitions of fairness, and their utility within different contexts is being actively debated [37, 41, 43, 48, 64]. Notable
formulations of fairness include individual fairness, which dictates that similar inputs must result in similar outputs, and
group fairness, which dictates that a particular subset of inputs must have a similar aggregate output to the whole. In this
paper, we view such notions of fairness as quantitative properties of decision-making programs.
Bias as a probabilistic property We think of decision-making algorithms as probabilistic programs, in the sense that
they are invoked on input drawn from a probability distribution, e.g., representing the demographics of some population.
Quantifying program bias becomes a matter of reasoning about probabilities of program executions.
Consider a hiring program P that takes as input a vector of arguments v representing a job applicant’s record. One of
the arguments vs in the vector v states whether the person is a member of a protected minority or not, and similarly vq in v
states whether the person is qualified or not. Evaluating P(v) returns a Boolean value indicating whether a person is hired.
Our goal may be to prove a group fairness property that is augmented with a notion of qualification—that the algorithm is
just as likely to hire a qualified minority applicant as it is for other qualified non-minority applicants. Formally, we state
this probabilistic condition as follows:
Pr[P(v) = true | vs = true ∧ vq = true]
>1−
Pr[P(v) = true | vs = false ∧ vq = true]
Here,  is a small prespecified value. In other words, the probability of hiring a person v, conditioned on them being
a qualified minority, is very close to (or greater than) the probability of hiring a person conditioned on them of being
1

2017/3/8

a qualified non-minority. The goal of this paper is to propose techniques for automatically proving whether a program
satisfies this kind of probabilistic properties.
Proving statements of the above form amounts to quantifying probabilities of return outcomes of the program. We
propose an automated verification technique that reduces the verification problem to that of computing the weighted
volume of the logical encoding of a program in real arithmetic. We then utilize a novel symbolic volume computation
algorithm that exploits the power of SMT solvers to integrate probability density functions over real regions defined by
program encodings. We show that our algorithm is guaranteed to converge to the exact values in the limit, thus resulting in
a sound and complete verification procedure. To our knowledge, this is one of the first probabilistic inference algorithms
for SMT with this expressivity and guarantees. We implement our algorithm in a tool called FairSquare, which we evaluate
on a number decision-making programs generated by a range of machine-learning algorithms from real-world data.
Contributions This paper makes a number of conceptual, algorithmic, and practical contributions:
– We present a verification technique for probabilistic programs that reduces the problem to a set of weighted volume
computation problems. We present a novel weighted-volume-computation algorithm, for formulas over real closed
fields, that utilizes an SMT solver as a black box, and we prove that it converges to the exact volume in the limit. To
our knowledge, this is one of the first probabilistic inference algorithms for SMT with this generality and guarantees.
(Sec. 4)
– We present an automated verification tool, FairSquare, and use it to quantify certain types of bias in a broad spectrum
of programs representing machine-learning classifiers generated from real-world datasets. Our evaluation demonstrates
the power of our technique and its ability to outperform state-of-the-art probabilistic program analyses. (Sec. 5)

2.

Overview and Illustration

Our problem setting is as follows: First, we are given a decision-making program Pdec . Second, we have a probabilistic
precondition defining a probability distribution over inputs of Pdec . We define the probability distribution operationally as
a probabilistic program Ppre , which we call the population model. Intuitively, the population model provides a probabilistic
picture of the population from which the inputs of Pdec are drawn. Third, we are given a quantitative postcondition post
that correlates the probabilities of various program outcomes. This postcondition can encode various properties relating to
program bias; intuitively, our goal is to prove the following triple:
{v ∼ Ppre }

r ← Pdec (v) {post}

Let us consider various possible instantiations of post. Feldman et al. [41] introduced the following definition, inspired
by Equality of Employment Opportunity Commission’s [6] recommendation in the US:
Pr[r = true | min(v) = true]
>1−
Pr[r = true | min(v) = false]
Assuming Pdec returns a Boolean value—indicating whether an applicant v is hired—this group fairness property states
that the selection rate from a minority group, min(v) = true, is as good as the selection rate from the rest of the population.
One can thus view this verification problem as proving a probabilistic property involving two sets of program traces: one
set where the input min(v) is true, and another where it is false. Alternatively, the above definition could be strengthened
with a lower bound on the ratio, so as to ensure that the selection rate of the two groups is similar (statistical parity).
Further, we could additionally condition on qualified applicants, e.g., if the job has some minimum qualification, we
do not want to characterize group fairness for arbitrary applicants, but only within the qualified subpopulation. Various
comparable notions of group fairness have been proposed and used in the literature, e.g., [33, 41, 72].
While the above definition is concerned with fairness at the level of subsets of the domain of the decision-making
program, individual fairness [37] is concerned with similar outcomes for similar elements of the domain. In our hiring
example, one potential formulation is as follows:
Pr[r1 = r2 | v 1 ∼ v 2 ] > 1 − 
In other words, we want to ensure that for any two individuals, if they are similar (∼), then we want them to receive similar
outcomes (r1 = r2 ) with a high probability. This is a hyperproperty—as it considers two copies of Pdec —and can be
encoded through self-composition [13]. This property is close in nature to differential privacy [36] and robustness [15, 24].
Of course, various definitions of fairness have their merits and their shortcomings, and there is an ongoing discussion
on this subject [7, 37, 41, 43, 48]. Our contribution is not to add to this debate, but to cast fairness as a quantitative property
of programs, and therefore enable automated reasoning about fairness of decision-making programs.
2

2017/3/8

A simple verification problem Consider the two programs in Figure 1(a). The program popModel is a probabilistic program
describing a simple model of the population. Here, a member of the population has three attributes, all of which are realvalued: (i) ethnicity; (ii) colRank, the ranking of the college the person attended (lower is better); and (iii) yExp, the years
of work experience a person has. We consider a person is a member of a protected group if ethnicity > 10; we call this
the sensitive condition. The population model can be viewed as a generative model of records of individuals—the more
likely a combination is to occur in the population, the more likely it will be generated. For instance, the years of experience
an individual has (line 4) follows a Gaussian (normal) distribution with mean 10 and standard deviation 5. Observe that
our model specifies that members of a protected minority will probably attend a lower-ranked college, as encoded in lines
5-6.
The program dec is a decision-making program that takes a job applicant’s college ranking and years of experience and
decides whether they get hired (the fairness target). The program implements a decision tree, perhaps one generated by a
machine-learning algorithm. A person is hired if they attended a top-5 college (colRank <= 5) or have lots of experience
compared to their college’s ranking (expRank > -5). Observe that dec does not access an applicant’s ethnicity.
Our goal is to establish whether the hiring algorithm dec discriminates against members of the protected minority.
Concretely, we attempt to prove the following property:
Pr[hire | min]
>1−
Pr[hire | ¬min]
where min is shorthand for the sensitive condition ethnicity > 10, and  is a small parameter set to 0.1 for illustration.
Despite the potential shortcomings of this group fairness property [37], its simple formulation serves well as an illustration
of our technique.
We can rewrite the above statement to eliminate conditional probabilities as follows:
Pr[hire ∧ min] · Pr[¬min]
>1−
Pr[hire ∧ ¬min] · Pr[min]

(1)

Therefore, to prove the above statement, we need to compute a value for each of the probability terms: Pr[hire ∧ min],
Pr[min], and Pr[hire ∧ ¬min]. (Note that Pr[¬min] = 1 − Pr[min].)
Notice that, to prove or disprove inequality 1, all we need are good enough bounds on the values of such probabilities
and not their exact values.
For the purposes of illustration, we shall focus our description on computing Pr[hire ∧ ¬min].
Probabilistic verification conditions To compute the probability Pr[hire ∧¬min], we need to reason about the composition
of the two programs, dec◦popModel. That is, we want to compute the probability that (i) popModel generates a non-minority
applicant, and (ii) dec hires that applicant. To do so, we begin by encoding both programs as formulas in the linear real
arithmetic theory of first-order logic. The process is analogous to that of standard verification condition (VC) generation
for loop-free program fragments.
First, we encode popModel as follows:
ϕpop ≡ ethnicity > 10 ⇒ colRank1 = colRank + 5
∧ ethnicity 6 10 ⇒ colRank1 = colRank

where subscripts are used to encode multiple occurrences of the same variable (i.e., SSA form). Note that assignments in
which values are drawn from probability distributions do not appear in the encoding—we shall address them later.
Second, we encode dec as follows (after simplification):
ϕdec ≡ expRank = yExpi − colRanki

∧ hire ⇐⇒ (colRanki 6 5 ∨ expRank > −5)

where variables with the superscript i are the input arguments to dec. Now, to encode the composition dec ◦ popModel, we
simply conjoin the two formulas—ϕpop and ϕdec —and add equalities between returns of popModel and arguments of dec.
ϕP ≡ ϕpop ∧ ϕdec ∧ yExpi = yExp ∧ colRanki = colRank1
Our goal is to compute the probability that a non-minority applicant gets hired. Formally, we are asking, what is the
probability that the following formula is satisfied?
ϕ ≡ ∃Vd . ϕP ∧ hire ∧ ethnicity 6 10
3

2017/3/8

(b)

1 define popModel()

Formula ' in R3
(blue faces are unbounded)

2

ethnicity ~ gauss(0,10)

3

colRank ~ gauss(25,10)

4

yExp ~ gauss(10,5)

5

if (ethnicity > 10)

6
7

}

The decision-making program takes an
applicant’s record and decides whether
to hire them. A person is hired if they
are from a top-5 college (colRank <= 5)
or have lots of experience compared to
their college rank (expRank > -5). Note
that this program does not access an
applicant’s ethnicity.

}

(a)
The population model defines a joint probability distribution on attributes of members of a population: (i) the rank of the
college a person attended (colRank), (ii) the
years of work experience they have (yExp),
and (iii) their ethnicity (ethnicity). Note
that colRank is influenced by a persons’s
ethnicity.

colRank

colRank + 5

return colRank, yExp

1 define dec(colRank, yExp)
2

expRank

3

if (colRank <= 5)

4
5
6
7
8
9

(c)

hire

yExp - colRank
true

elif (expRank > -5)
hire

true

else
hire

false

return hire

Underapproximation of '
as a union of hyperrectangles

colRank

ethnicity

FairSquare ratio computation
on dec and popModel (unfair)

FairSquare ratio computation
on modified dec (fair)

ratio upper/lower bound

(e)

ratio upper/lower bound

(d)

yExp

# of iterations of algorithm

# of iterations of algorithm

Figure 1. Simple illustrative example

Vd is the set of variables that are not probabilistically assigned to, that is, all variables other than the three variables
Vp = {ethnicity, colRank, yExp}. Intuitively, by projecting out all non-probabilistic variables, we get a formula ϕ whose
models are the set of all probabilistic samplings that lead to a non-minority applicant being generated and hired.
Weighted volume computation To compute the probability that ϕ is satisfied, we begin by noting that ϕ is, geometrically,
a region in R3 , because it has three free real-valued variables, Vp . The region ϕ is partially illustrated in Figure 1(b).
Informally, the probability of satisfying ϕ is the probability of drawing values for the dimensions in Vp that end up falling
in the region ϕ. Therefore, the probability of satisfying ϕ is its volume in R3 , weighted by the probability density of each
of the three variables. Formally:
R
Pr[hire ∧ ¬min] = ϕ pe py pc dVp
where, e.g., pe is the probability density function of the distribution gauss(0,10)—the distribution from which the
value of ethnicity is drawn in line 2 of popModel. Specifically, pe is a function of ethnicity, namely, pe (ethnicity) =
2

ethnicity
1
√
e− 200
10 2π

.
The primary difficulty here is that the region of integration is specified by an arbitrary SMT formula over an arithmetic
theory. So, how do we compute a numerical value for this integral? We make two interdependent observations: (i) if
4

2017/3/8

the formula represents a hyperrectangular region in Rn —i.e., a box—then integration is typically simple, due to the
constant upper/lower bounds of all dimensions; (ii) we can symbolically decompose an SMT formula into an (infinite) set
of hyperrectangles.
Specifically, given our formula ϕ, we construct a new formula, ϕ , where each model m |= ϕ corresponds to a
hyperrectangle that underapproximates ϕ. Therefore, by systematically finding disjoint hyperrectangles inside of ϕ and
computing their weighted volume, we iteratively improve a lower bound on the exact weighted volume of ϕ. Figure 1(c)
shows a possible underapproximation of ϕ composed of four hyperrectangles. Sec. 4 formalizes this technique and proves
its convergence for decidable arithmetic theories.
Proofs of group fairness We demonstrated how our technique reduces the problem of computing probabilities to weighted
volume computation. Figure 1(d) illustrates a run of our tool, FairSquare, on this example. FairSquare iteratively improves
lower and upper bounds for the probabilities in the ratio, and, therefore, the ratio itself. Observe how the upper bound (red)
of the ratio is decreasing and its lower bound (blue) is increasing. This example is not group fair for  = 0.1, since the
upper bound goes below 0.9.
Recall that applicants of a protected minority tend to attend lower-ranked colleges, as defined by popModel. Looking
at dec, we can point out that the cause for unfairness is the importance of college ranking for hiring. Let us attempt to fix
this by modifying line 2 of dec to expRank ← 5*yExp - colRank. In other words, we have made the hiring algorithm value
an applicant’s job experience way more than college ranking. The run of FairSquare on the modified dec is illustrated in
Figure 1(e), where the lower bound on the ratio exceeds 0.9, thus proving our fairness property.

3.

Probabilistic Programs and Verification

We formally define programs and present a general framework for stating and verifying probabilistic properties.
3.1

Program model and semantics

Programs A program P is a sequence of statements S:
S := V ← E
|V ∼D
| if B then S else S
| SS

assignment statement
probabilistic assignment
conditional
sequence of statements

where V is the set of real-valued variables that can appear in P, e ∈ E is an arithmetic expression over variables in V ,
and b ∈ B is a Boolean expression over variables in V . A probabilistic assignment is made by sampling from a probability
distribution p ∈ D. A probability distribution can be, for example, a Gaussian distribution, denoted by gauss(µ, σ), where
µ, σ ∈ R are the mean and standard deviation of the Gaussian. We shall restrict distributions to be univariate and with
constant parameters, e.g., mean and standard deviation of a Laplacian or Gaussian—that is, we assume independence of
probabilistic assignments. Given a probabilistic assignment x ∼ p, we shall treat p(x) as a probability density function
(PDF) of the distribution from which the value assigned to x is drawn. For instance, if the distribution is gauss(0,1), then
x2
p(x) = √12π e− 2 .
We use v i to denote a vector of input variables of P, and v o to denote a vector of output variables of P; these variables
appear in V and denote the arguments and returns of P. We say that a program is closed if it has no inputs, i.e., v i is
empty. We shall refer to the following subsets of V .
– Vp ⊆ V is the set of probabilistic variables: those that get assigned to in probabilistic assignments.

– Vd = V \ Vp is the set of deterministic variables: those that do not appear in probabilistic assignments.
This simple language can be used to describe decision-making programs typical machine-learning classifiers such as
decision trees, support vector machines, Bayesian networks, neural networks, as well as loop-free probabilistic programs
(loops with constant bounds can be unrolled).1 As demonstrated in Sec. 2, the same language is used to define probabilistic
preconditions programmatically.
Operational semantics Typically, the state s : V → R of the program is defined as a valuation function from variables
in V to values in R. In a probabilistic setting, however, we need to maintain an additional state that dictates values drawn
from probability distributions. Following standard semantics of probabilistic programs [52], we assume a finite sequence
of independent random variables. The semantics of an execution is thus defined for a fixed sequence of values ω of these
variables. Informally, it is as if we performed all sampling before the program executes and stored the results in a sequence
for use whenever we encounter a probabilistic assignment. For the full semantics, refer to Appendix A.
1 Unbounded

loops could be handled through iterative unrolling as in [66].

5

2017/3/8

VC - ASN

hx = JeK, ∅i  x ← e

hϕ1 , D1 i  S1

VC - PASN

htrue, {pi }i  xi ∼ pi

hϕ2 , D2 i  S2

VC - SEQ

hϕ1 ∧ ϕ2 , D1 ∪ D2 i  S1 S2

hϕ1 , D1 i  S1

hϕ2 , D2 i  S2

VC - COND

hite(JbK, ϕ1 , ϕ2 ), D1 ∪ D2 i  if b then S1 else S2

Figure 2. Probabilistic verification condition generation. Above, ite(a, b, c) , (a ⇒ b) ∧ (¬a ⇒ c).

3.2

Programs and volume computation

Following Chistikov et al. [27], we reduce the problem of computing the probability that the program terminates in a state
satisfying ϕ to weighted volume computation (WVC) over formulas describing regions in Rn . In what follows, we begin
by formalizing the WVC problem.
Volume of a formula We will use LT to denote first-order formulas in some real arithmetic theory T . Specifically, we
consider two decidable theories: linear real arithmetic and the strictly richer real closed fields—Boolean combinations of
polynomial inequalities. Given a formula ϕ ∈ LT , a model m of ϕ, denoted by m |= ϕ, is a point in Rn , where n is the
number of free variables of ϕ. Thus, we view ϕ as a region in Rn , i.e., ϕ ⊆ Rn . We use Xϕ = {x1 , . . . , xn } to denote the
free variables of ϕ.
R
The (unweighted) volume of a formula ϕ is ϕ 1 dXϕ where dXϕ is short for dx1 dx2 . . . dxn . For example, if ϕ is in
R
R2 , then ϕ 1 dXϕ is the area of ϕ.
Weighted volume of a formula We now define the weighted volume of a formula. We assume we are given a pair (ϕ, D),
where ϕ ∈ LT and D = {p1 , . . . , pn } is a set of probability density functions such that each variable xi ∈ Xϕ is associated
with a density function pi (xi ) of the probability distribution of its values. The weighted volume of ϕ with respect to D,
denoted by VOL(ϕ, D), is defined as follows:
R Q
xi ∈Xϕ pi (xi ) dXϕ
ϕ
Example 1. Consider the formula ϕ ≡ x1 + x2 > 0, and let D = {p1 , p2 }, where p1 and p2 are the PDF of the Gaussian
distribution with mean 0 and standard deviation 1. Then,
R
VOL(ϕ, D) = x +x >0 p1 (x1 )p2 (x2 ) dx1 dx2 = 0.5
1
2
Intuitively, if we are to randomly draw two values for x1 and x2 from the Gaussian distribution, we will land in the region
x1 + x2 > 0 with probability 0.5.

Probabilistic verification conditions Recall that our goal is to compute the probability of some predicate ϕ at the end of a
program execution, denoted Pr[ϕ]. We now show how to encode this problem as weighted volume computation. First, we
encode program executions as a formula ϕP . The process is similar to standard verification condition generation (as used
by verification [11] and bounded model checking tools [30]), with the difference that probabilistic assignments populate a
set D of probability density functions.
Figure 2 inductively defines the construction of a probabilistic verification condition for a program P, denoted by a
function PVC(P), which returns a pair hϕP , Di. Without loss of generality, to simplify our exposition, we assume programs
are in static single assignment (SSA) form [31]. Given a Boolean expression b, the denotation JbK is the same expression
interpreted as an LT formula. The same applies to arithmetic expressions e. For example, Jx + y > 0K , x + y > 0.
Intuitively, the construction generates a formula ϕP that encodes program executions, treating probabilistic assignments
as non-deterministic, and a set D of the PDFs of distributions in probabilistic assignments (rule VC - PASN).
Now, suppose we are given a closed program P and a Boolean formula ϕ over its output variables. Then,
Pr[ϕ] = VOL(∃Vd . ϕP ∧ ϕ, D)
That is, we project out all non-probabilistic variables from ϕP ∧ ϕ and compute the weighted volume with respect to the
densities pi ∈ D. Intuitively, each model m of ∃Vd . ϕP ∧ ϕ corresponds to a sequence of values drawn in probabilistic
assignments in an execution of P. We note that our construction is closely related to that of Chistikov et al. [27], to which
we refer the reader for a measure-theoretic formalization.
6

2017/3/8

1:
2:
3:
4:
5:
6:
7:
8:
9:

function VERIFY(Ppre , Pdec , post)
hϕpre , Dpre i ← PVC(Ppre )
hϕdec , Ddec i ← PVC(Pdec )
pre
hϕP , Di ← hϕpre ∧ ϕdec ∧ v dec
i = v o , Dpre ∪ Ddec i
pre
dec
Vd ← Vd ∪ Vd
m←∅
for Pr[ϕ] ∈ post do
m ← m[Pr[ϕ] 7→ VOL(∃Vd . ϕP ∧ ϕ, D)]
return m |= post
Figure 3. Abstract verification algorithm

Example 2. Consider the following closed program P
x ~ gauss(0,2); y ~ gauss(-1,1); z

←

x + y

where z is the return variable. Using the encoding in Figure 2, we compute the pair hϕP , Di  P, where ϕP , z = x + y
and D = {px , py }, where px and py are the PDFs of the two distributions from which values of x and y are drawn.
Suppose that we would like to compute the probability that z is positive when the program terminates: Pr[z > 0]. Then,
we can compute the following weighted volume: VOL(∃z. ϕP ∧ z > 0, D), which is equal to 0.32736.

3.3

Probabilistic verification problems

We now define probabilistic verification problems and present an abstract verification algorithm that assumes the existence
of an oracle for weighted volume computation.
Verification problems A verification problem is a triple (Ppre , Pdec , post), where
– Ppre , called the probabilistic precondition, is a closed program over variables V pre and output variables v pre
o .

– Pdec , called the decision-making program, is an open program over variables V dec ; its input arguments are v dec
i , with
pre
dec
dec
pre
.
(We
assume
that
V
∩
V
=
∅.)
|;
and
its
output
variables
are
v
|
=
|v
|v dec
o
o
i

– post is a probabilistic postcondition, which is a Boolean expression over probabilities of program outcomes. Specifically, post is defined as follows:
post ∈ PExp := Prob > Prob | PExp ∨ PExp
| PExp ∧ PExp | ¬PExp

Prob := Pr[ϕ] | Pr[ϕ | ϕ] | c | Prob

Prob

∈ {+, −, ÷, ×} c ∈ R

where ϕ ∈ LT are formulas over input and output variables of Pdec . For example, post might be of the form
Pr[x > 0] > 0.5 ∧ Pr[y + z > 7] > Pr[t > 5]
The goal of verification is to prove that post is true for the program Pdec ◦ Ppre , i.e., the composition of the two programs
where we first run Ppre to generate an input for Pdec . Since Ppre is closed, the program Pdec ◦ Ppre is also closed.
Verification algorithm We now describe an idealized verification algorithm that assumes the existence of an oracle for
weighted volume computation. The algorithm, VERIFY, shown in Figure 3, takes a verification problem and returns whether
the probabilistic postcondition holds.
VERIFY begins by encoding the composition of the two programs, Pdec ◦Ppre , as the pair hϕP , Di and adds the constraint
pre
v dec
i = v o to connect the outputs of Ppre to the inputs of Pdec (recall the example from Sec. 2 for an illustration). For each
term of the form Pr[ϕ] appearing in post, the algorithm computes its numerical value and maintains it in a map m. If m
satisfies the post—i.e., by replacing all terms Pr[ϕ] with their values in m—then the postcondition holds.

4.

Symbolic Weighted Volume Computation

In this section, we describe our weighted volume computation algorithm. Recall that, Rgiven
Q a formula ϕ and a set D
defining the PDFs of the distributions of free variables, our goal is to evaluate the integral ϕ xi ∈Xϕ pi (xi ) dXϕ .
Existing techniques In general, there is no systematic technique for computing an exact value for such an integral.
Moreover, even simpler linear versions of the volume computation problem, not involving probability distributions,
are #P-hard [39]. Existing techniques suffer from one or more of the following: they (i) restrict ϕ to a conjunction of
7

2017/3/8

vol ← 0

Ψ ← ϕ

HDECOMP

m |= Ψ
m

vol ← vol + VOL(H , D)
where block (H m ) ≡

Ψ ← Ψ ∧ block (H m )

_
x∈Xϕ

HSAMPLE

ux < Hlm (x) ∨ lx > Hum (x)

Figure 4. SYMVOL: weighted volume computation algorithm

linear inequalities [34, 66], (ii) restrict integrands to polynomials or simple distributions [16, 17, 28, 34], (iii) compute
approximate solutions with probabilistic guarantees [17, 27, 28, 70], (iv) restrict ϕ to bounded regions of Rn [27, 28], or
(v) have no convergence guarantees, e.g., computer algebra tools that find closed-form solutions [2, 3, 44]. (See Sec. 6 for
details.)
Symbolic weighted volume computation Our approach is novel in its generality and its algorithmic core. The following
are the high-level properties of our algorithm:
1. Our approach accepts formulas in the decidable yet rich theory of real closed fields: Boolean combinations of
polynomial inequalities.
2. Our approach imposes no restrictions on the form of the PDFs, only that we can evaluate the cumulative distribution
functions2 (CDFs) associated with the PDFs in D.

3. Our approach is guaranteed to converge to the exact value of the weighted volume in the limit, allowing us to produce
a sound and complete verification procedure.
At the algorithmic level, our approach makes the following contributions:
1. Our approach exploits the power of SMT solvers and uses them as a black box, allowing it to directly benefit from future
advances in solver technology.
2. Our approach employs the idea of dividing the space into rectangular regions that are easy to integrate over. While this
age-old idea has been employed in various guises in verification [9, 21, 53, 66], we utilize it in a new symbolic way to
enable volume computation over SMT formulas.
3. Our approach introduces a novel technique for approximately encoding PDFs as formulas, and using them to guide the
SMT solver towards making large leaps to the exact solution.
4.1

Weighted Volume Computation Algorithm

To compute the integral over the region ϕ, we exploit the observation that if ϕ is a hyperrectangular region, i.e., an ndimensional rectangle in Rn , then we can evaluate the integral, because each dimension has constant lower and upper
bounds. For instance, consider the following formula representing a rectangle in R2 :
ϕ ≡ 0 6 x1 6 100 ∧ 4 6 x2 6 10
R

p (x )p (x ) dx1 dx2
ϕ 1 1 2 2
R 100
R 10
=( 0 p1 (x1 ) dx1 )( 4 p2 (x2 ) dx2 )

The following holds:

=(F1 (10) − F1 (4))(F2 (100) − F2 (0))
Rx
where Fi is the CDF of pi (xi ), which is Fi (x) = −∞ pi (t) dt. That is, we independently compute the integral along each
dimension of the rectangle, and take the product. This holds since we assume all variables are independently sampled.
Our algorithm is primarily composed of two steps: First, the hyperrectangular decomposition phase represents the
formula ϕ as a set of hyperrectangles. Note that this set is likely to be infinite. Thus, we present a technique for defining
all hyperrectangles that lie in ϕ symbolically as a formula ϕ , where each model of ϕ corresponds to a hyperrectangle
2 The

cumulative distribution function of a real-valued random variable X is the function f : R → R, such that f (x) = Pr[X 6 x].

8

2017/3/8

high probability

Infinitely many
rectangles are required to
cover this triangular
polytope

Point rectangles
have no volume

(a)

Infinitely many
rectangles in low
probability regions
have negligible
volume

high probability

(b)

Figure 5. R2 view of (a) hyperrectangular decomposition and (b) hyperrectangle sampling, where density is concentrated in the top-left corner

that lies inside the region ϕ. Second, after characterizing the set ϕ of all hyperrectangles in ϕ, we can iteratively sample
hyperrectangles in ϕ, which can be done using an off-the-shelf SMT solver to find models of ϕ . For each hyperrectangle
we sample, we compute its weighted volume and add it to our current solution. Therefore, the current solution maintained
by the algorithm is the weighted volume of an underapproximation of ϕ—that is, a lower bound on the exact weighted
volume of ϕ.
Hyperrectangular decomposition We begin by defining hyperrectangles as special formulas.
Definition 1V(Hyperrectangles and their weighted volume). A formula H ∈ LT is a hyperrectangle if it can be written
in the form x∈XH cx 6 x 6 c0x where cx , c0x ∈ R are the lower and upper bounds of dimension x. We use Hl (x) and
Hu (x) to denote the lower and upper bounds of x in H.
The weighted volume of H, given a set D, is as follows:
VOL(H, D)

=

Q

xi ∈XH

R Hu (xi )
Hl (xi )

pi (xi ) dxi

W
Ideally, we would take a formula ϕ and rewrite it as a disjunction of hyperrectangles H, but this disjunction is
most likely an infinite one. To see why this is the case, consider the simple formula representing a triangular polytope in
Figure 5(a). Here, there is no finite number of rectangles whose union is the full region in R2 enclosed by the triangle.
While the number of hyperrectangles enclosed in ϕ is infinite, we can characterize them symbolically using universal
quantifiers, as shown by Li et al. [53]. Specifically, we define the hyperrectangular decomposition of ϕ as follows:
Definition 2 (Hyperrectangular decomposition). Given a formula ϕ, its hyperrectangular decomposition ϕ is:
 ^

 ^

ϕ ≡
lx 6 ux ∧ ∀Xϕ .
lx 6 x 6 u x ⇒ ϕ
x∈Xϕ

x∈Xϕ

where lx , ux are fresh free variables introduced for each x ∈ Xϕ , and ∀Xϕ is short for ∀x1 , . . . , xn , for xi ∈ Xϕ .
Given a model m |= ϕ , we say that H m is the hyperrectangle induced by m, as defined below:
^
m(lx ) 6 x 6 m(ux )
Hm ≡

Intuitively, ϕ characterizes every possible hyperrectangle that is subsumed by ϕ. The idea is that the hyperrectangle
H m induced by each model m of ϕ is subsumed by ϕ, that is, H m ⇒ ϕ. The following example illustrates this process.
x∈Xϕ

Example 3. Consider the formula ϕ ≡ x > y ∧ y > 0, illustrated in Figure 6 as a gray, unbounded polyhedron. The
formula ϕ , after eliminating the universal quantifier, is:
lx 6 ux ∧ ly 6 uy ∧ ly > 0 ∧ lx > uy

Figure 6 shows two models m1 , m2 |= ϕ and their graphical representation as rectangles H m1 , H m2 in R2 . Observe
that both rectangles are subsumed by ϕ.

Hyperrectangle sampling Our symbolic weighted volume computation algorithm, SYMVOL, is shown in Figure 4 as two
transition rules. Given a pair (ϕ, D), the algorithm maintains a state consisting of two variables: (i) vol , the current lower
bound of the weighted volume, and (ii) Ψ, a constraint that encodes the remaining rectangles in the hyperrectangular
decomposition of ϕ.
9

2017/3/8

y

Two examples of models
inducing rectangles

m1 |= ⌧'
lx = 8
ux = 10
ly = 5
uy = 7

m2 |= ⌧'
lx = 5
ux = 7
ly = 0
uy = 5

'⌘x>y^y >0

7

H m1

5

H m2
5

Figure 6. Illustration of models of ϕ
(a)

(b)

7 8

10 x

(c)

Figure 7. Three ADFs (gray) of a Gaussian PDF (red) with mean 0 and a standard deviation 3. ADFs are: (a) fine-grained; (b) coarse; (c) uniform.

The algorithm is presented as guarded rules. Initially, using the rule HDECOMP, vol is set to 0 and Ψ is set to ϕ . The
algorithm then proceeds by iteratively applying the rule HSAMPLE. Informally, the rule HSAMPLE is used to find arbitrary
hyperrectangles in ϕ and compute their weighted volume. Specifically, HSAMPLE finds a model m of Ψ, computes the
weighted volume of the hyperrectangle H m induced by m, and adds the result to vol .
To maintain soundness, HSAMPLE ensures that it never samples two overlapping hyperrectangles, as otherwise we
would overapproximate the volume. To do so, every time a hyperrectangle H m is sampled, we conjoin an additional
0
constraint to Ψ—denoted block (H m ) and defined in Figure 4—that ensures that for all models m0 |= Ψ, H m does not
0
overlap with H m , i.e., H m ∧ H m is unsatisfiable. Informally, the block (H m ) constraint specifies that any newly sampled
hyperrectangle should be to the left or right of H m for at least one of the dimensions.
Lower and upper bounds The following theorem states the soundness of SYMVOL: it maintains a lower bound on the
exact weighted volume.
Theorem 1 (Soundness of SYMVOL). The following is an invariant of SYMVOL(ϕ, D): vol 6 VOL(ϕ, D).
It follows from the above theorem that we can use SYMVOL to compute an upper bound on the exact volume.
Specifically, because we are integrating over PDFs, we know that VOL(ϕ, D) + VOL(¬ϕ, D) = 1. Therefore, by using
SYMVOL to compute the weighted volume of ¬ϕ, we get an upper bound on the exact volume of ϕ.

Corollary 1 (Computing upper bounds). The following is an invariant of SYMVOL(¬ϕ, D): 1 − vol > VOL(ϕ, D)
4.2

Density-directed Sampling

While the SYMVOL algorithm is sound, it provides no progress guarantees. Consider, for example, a run that only samples
unit hyperrectangles, i.e., points in Rn ; the volume of a point is 0, therefore, we will never compute any volume.
Alternatively, the algorithm might diverge by sampling hyperrectangles in ϕ that appear in very low probability density
regions. These two scenarios are illustrated in Figure 5(b) on a triangular polytope in R2 .
Ideally, the rule HSAMPLE would always find a model m yielding the hyperrectangle H m with the largest weighted
volume. Finding such a model amounts to solving the optimization problem:
R H m (x )
Q
arg max xi ∈Xϕ H mu (xii) pi (xi ) dxi
l

m|=Ψ

From a practical perspective, there are no known tools or techniques for finding models of first-order formulas that
maximize such complex objective functions, with integrals over arbitrary probability density functions.
However, we make theR key observation that if p(x) is a step function—i.e., piecewise constant—then we can symbolically encode the integral p(x) dx in linear arithmetic.
R As such, we propose to (i) approximate each density function p(x)
with a step function step(x), (ii) encode the integrals step(x) dx as linear arithmetic formulas, and (iii) direct sampling
towards hyperrectangles that maximize these integrals, thus finding hyperrectangles of large volume.
10

2017/3/8

vol ← 0
m |= Ψ ∧

^

Ψ ← ϕ

HDECOMP

lb ← 1

∃δxi . step φ
i (xi ) ∧ δxi > lb

xi ∈Xϕ

vol ← vol + VOL(H m , D)

Ψ ← Ψ ∧ block (H m )

HSAMPLE

DECAY

lb ← λ ∗ lb
Figure 8. ADF - SYMVOL: ADF-directed volume computation

Approximate density functions and linear encodings We begin by defining approximate density functions (ADFs).
Definition 3 (Approximate density functions). An approximate density function step(x) is of the following form:
(
ci , x ∈ [ai , bi ) for 1 6 i 6 n
step(x) =
0, otherwise
where ci , ai , bi ∈ R, ci > 0, and all [ai , bi ) are disjoint.

We now show how to encode a formula step φ (x) over the free variables δx , lx , ux , where for any model m |= step φ (x),
R m(u )
the value m(δx ) is the area under step(x) between m(lx ) and m(ux ), i.e.: m(δx ) = m(lxx) step(x) dx. Intuitively, the
value of this integral is the sum of the areas of each bar in step(x), restricted to[m(lx ), m(ux )].
Definition 4 (Encoding area under an ADF). Given an ADF step(x), we define step φ (x) as follows:
Pn
step φ (x) ≡ δx = i=1 ci · [ai , bi ) ∩ [lx , ux ]
The constraint step φ (x) is directly expressible in linear arithmetic, since
[ai , bi ) ∩ [lx , ux ] = max(min(bi , ux ) − max(ai , lx ), 0)
The finite sum in step φ (x) computes the size of the intersection of [lx , ux ] with each interval [ai , bi ) in step(x), and
multiplies the intersection with ci , the value of the step in that interval.
ADF -directed Volume Computation We now present the algorithm ADF - SYMVOL (Figure 8), an extension of our volume
computation algorithm SYMVOL that uses ADFs to steer the sampling process. The ADFs are only used for guiding the
rule HSAMPLE towards dense hyperrectangles, and thus do not affect soundness of the volume computation. For example,
Figure 7 shows three approximations of a Gaussian; all three are valid approximations. In Section 5, we discuss the impact
of different ADFs on performance.
Formally, we create a set of ADFs A = {step 1 , . . . , step n }, where, for each variable xi ∈ Xϕ , we associate the ADF
step i (xi ). The rule HSAMPLE now encodes step φi (xi ) and attempts to find a hyperrectangle such that for each dimension
x, δx is greater than some lower bound lb, which is initialized to 1. Of course, we need to keep reducing the value lb as
we run out of hyperrectangles of a given volume. Therefore, the rule DECAY is used to shrink lb using a fixed decay rate
λ ∈ (0, 1) and can be applied when HSAMPLE fails to find a sufficiently
large hyperrectangle.
Q
Note that, ideally, we would look for a model m such that x∈Xϕ δx is maximized, thus, finding the hyperrectangle
with the largest weighted volume with respect to the ADFs. However, this constraint is non-linear. To lower the complexity
of the problem to that of linear arithmetic, we set a decaying lower bound and attempt to find a model where each δx is
greater than the lower bound.
4.3

Convergence of ADF - SYMVOL

We now discuss the convergence properties of ADF - SYMVOL. Suppose we are given a formula ϕ, a set D, and a set A. Let
R ⊂ Rn be the region where all the ADFs in A are non-zero. We will show that
- SYMVOL converges, in the limit, to
R ADFQ
the exact weighted volume restricted to R; that is, ADF - SYMVOL converges to ϕ∩R xi ∈Xϕ pi (xi ) dXϕ .
Note the fascinating part here is that we do not impose any restrictions on the ADFs: they do not have to have any
correspondence with the PDFs they approximate; they need only be step functions. Of course, in practice, the quality of the
approximation dictates the rate of convergence, but we delay this discussion to Section 5.
The following theorem states convergence of ADF - SYMVOL; it assumes that HSAMPLE is applied iteratively and DECAY
is only applied when HSAMPLE cannot find a model.
11

2017/3/8

Theorem 2 (Convergence to R). Assume ADF - SYMVOL is run on (ϕ, D) and a set of ADFs
R A that
Q are non-zero for
R ⊂ Rn . Let vol i be the value of vol after i applications of HSAMPLE. Then, limi→∞ vol i = ϕ∩R xi ∈Xϕ pi (xi ) dXϕ .
Note that the above theorem directly gives us a way to approach the exact volume. Specifically, by performing runs
of ADF - SYMVOL on subsets in an infinite partition of Rn induced by the ADFs, we can ensure that the sum over the
ADF - SYMVOL processes approaches the exact volume in the limit. For all i, let Ai be a set of ADFs corresponding to an
ADF - SYMVOL process Pi , where Ri ⊂ Rn is the non-zero region of Ai . We require an infinite set of Pi to partition Rn :
S∞
(i) for all i 6= j, Ri ∩ Rj = ∅, and (ii) i=1 Ri = Rn . The following theorem formalizes the argument:

Theorem 3 (Convergence). Let P1 , P2 , . . . be ADF - SYMVOL processes that partition Rn . Assume a fair serialization
where each Pi performs HSAMPLE infinitely often, and let vol n be the total computed volume across all Pi after n
successful calls to HSAMPLE. Then, limn→∞ vol n = VOL(ϕ, D).

5.

Implementation and Evaluation

In this section, we present a case study on the bias of decision-making programs. In particular, we verify the group fairness
property with qualification, which requires that a program P satisfies the condition
Pr[P(v) = true | min(v) = true ∧ qual(v) = true]
> 0.85
Pr[P(v) = true | min(v) = false ∧ qual(v) = true]
We implement our algorithms in a tool called FairSquare and evaluate its ability to prove or disprove the property
for an assortment of programs; furthermore, we evaluate the effect of various parameters on the performance of the
implementation, and we compare the applicability of FairSquare to other probabilistic inference tools.
5.1

Implementation

We implemented our presented algorithms in a new tool called FairSquare, which employs Z3 [35] for SMT solving and
Redlog [4] for quantifier elimination. To compute the group fairness ratio, we decompose the conditional probabilities
into four joint probabilities: Pr[P(v) ∧ min(v) ∧ qual(v)], Pr[P(v) ∧ ¬min(v) ∧ qual(v)], Pr[min(v) ∧ qual(v)], and
Pr[¬min(v) ∧ qual(v)]. FairSquare computes lower bounds for eight quantities: the weighted volume of each of these
probabilities and also their negations, since a lower bound on the negation provides an upper bound on the positive form.
A round of sampling involves (i) obtaining a sample (hyperrectangle) for each of these eight quantities, (ii) computing
these samples’ weighted volumes, (iii) updating the bounds on the group fairness ratio, and (iv) checking if bounds are
precise enough to conclude fairness or unfairness. Rounds of sampling are performed until a proof is found or a timeout
is reached.
Sample maximization A key optimization is the maximization of hyperrectangles obtained during sampling. We use
Z3’s optimization capability to maximize and minimize the finite bounds of all hyperrectangles, while still satisfying the
formula Ψ (in Figures 4 and 8). This process is performed greedily by extending a hyperrectangle in one dimension at a
time to find a maximal hyperrectangle. If a dimension extends to infinity, then we drop that bound, thus resulting in an
unbounded hyperrectangle.
5.2

Benchmarks

We trained a variety of machine-learning models on a popular income dataset [1] used in related research on algorithmic
fairness [22, 41, 72] to predict whether a person has a low or high income; suppose, for example, these programs would be
used to determine the salary of a new employee: high (> $50,000) or low. We would like to verify whether salary decisions
are fair to qualified female employees. Using the Weka machine learning suite [5], we learned 11 different decision-making
programs (see, e.g., Bishop’s textbook [20] for background) which are listed in Figure 9:
– Four decision trees, named DTn , where n is the number of conditionals in the program. The number of variables and the
depth of the tree each varies from 2 to 3.
– Four support vector machines with linear kernels, named SVMn , where n is the number of variables in the linear
separator.
– Three neural networks using rectified linear units [60], named NNn,m , where n is the number of input variables, and m
is the number of nodes in the single hidden layer.
As we will show in the next section, some of these programs do not satisfy group fairness. We introduced modifications
α
of DT16 and SVM4 , called DTα
16 and SVM 4 , that implement rudimentary forms of affirmative action for female applicants.
For DTα
,
there
is
a
15%
chance
it
will
flip
a decision to give the low salary; for SVMα
16
4 , the linear separator is moved to
increase the likelihood of hiring.
12

2017/3/8

7
3
7
3
7
7
7
3
7

0.02
35.3

3
3
0.03
674.7

634 TO
78 21.6
62 27.8

2.4
0.8
2.0

442

10.0

TO

0.09
3.03

3
3
0.00
5.24

#

Vol

QE

Res

#

Vol

QE

12
2.2
0.9
7
38 52.3 11.4 3
22 15.3 6.3
7
34 32.0 8.2 3
113 178.9 94.3 7
10
3.7
1.7
7
18 13.3 3.1
7
22 15.7 3.2 3
10 12.2 6.3 TOq

18
6.6
2.2
73 130.9 33.6
22 38.2 14.3
40 91.0 19.4
406 484.0 222.4
10 10.8
6.2
14 33.7 20.1
14 33.4 63.2
TO

434
TO
12.8 TOq
466 456.1 3.4 3
238 236.5 7.2 3

154 132.9
174 233.5

34

TO

55.9 TOq

-

-

TO

7.2
18.2
TO

NN

10
8
6
4
2
0

FairSquare

PSI

[44]

VC

BN2

Res

0.5
1.4
2.0
3.0
9.8
0.6
0.8
0.8
1.3

Ind

0.74

QE

1.3
4.2
7.7
5.1
63.5
2.6
2.7
3.0
8.5

SVM

BN1

NN 3,2

Vol

10
20
21
18
55
10
10
10
10

DT

BN2

0.79
0.65
0.67

SVM 6

#

3
3
3
3
3
3
3
3
3

Ind

NN 2,1
NN 2,2

Res

BN1

0.79
0.71
0.79
0.76
0.82
0.79
0.79
0.78
0.79

12

Bayes Net 2

BN2

DT4
DT14
DT16
α
DT16
DT44
SVM 3
SVM 4
α
SVM 4
SVM 5

Population Model
Bayes Net 1

Independent

Ind

Acc

BN1

Decision
program

[66]

Figure 9. (Left) Results of FairSquare applied to 39 group fairness problems. Res indicates 3 for fair; 7 for unfair. Vol indicates the time (s) of the
sampling procedure; # is the number of SMT calls. If sampling timed out (900s), Res denotes the latest bounds on the fairness ratio. QE indicates the time
(s) of the quantifier elimination procedure used prior to sampling; if this timed out (900s), we could not perform any sampling, denoted by a TOq for Res.
Acc is the training set accuracy for each of the programs. (Right) Comparison of the number of benchmarks that FairSquare, PSI [44], and VC [66] were
able to solve.

Additionally, we used three different probabilistic population models, programs that define the probabilistic inputs, that
were inferred from the same dataset: (i) a set of independently distributed variables, (ii) a Bayesian network using a simple
graph structure, and (iii) the same Bayesian network, but with an integrity constraint in the form of an inequality between
two of the variables. Note that the first model is a trivial case: since there is there is no dependence between variables, all
programs will be fair; this simplicity serves well as a baseline for our evaluation. The Bayesian models permit correlations
between the variables, allowing for more subtle sources of fairness or unfairness. The benchmarks we use are derived from
each combination of population models with decision-making programs (see Appendix C for an example).
5.3

Evaluation

In this section, we discuss the ability of FairSquare to verify qualified group fairness or unfairness for 39 problems, as
summarized in Figure 9. In these problems, min(v) = true when the applicant is female, and qual(v) = true in two
different scenarios: first, we consider the case when qual is tautologically true, and second, when the applicant is at least
18 years of age. Figure 9 shows only the former case, as while the numbers are different for the two cases, the qualitative
results are quite similar. See Appendix D for the full table of results. We fix  = 0.15. To guide volume computation, all
Gaussian distributions with mean µ and variance σ 2 use ADFs with 5 equal-width steps spanning (µ − 3σ 2 , µ + 3σ 2 )—
analogous to Figure 7(a). FairSquare was able to solve 32 of the 39 problems, proving 21 fair and 11 unfair, as shown in
the left table of Figure 9.
Consider the results for DT4 : FairSquare proved it fair with respect to the independent population model after 0.5
seconds of an initial quantifier elimination procedure and 1.3 seconds of the actual volume computation algorithm, which
required 14 SMT queries. The more sophisticated Bayesian network models took longer for sampling, but due to the
correlations between variables, were proved unfair.
In contrast, consider the results for DT44 under the Bayes Net 1 population model: FairSquare was unable to conclude
fairness or unfairness after 900 seconds of volume computation (denoted by TO in the Vol column). The lower and upper
bounds of the fairness ratio it had computed at that time are listed in the Res column: in this case, the value of the fairness
ratio is within [0.70, 0.88], which is not precise enough for the  = 0.15 requirement (but would be precise enough for 
outside of [0.12, 0.30]).
In general, all conclusive results using the independent population model were proved to be fair, as expected, but many
are unfair with respect to the clusters and Bayes net models because of the correlations they capture. This difference
illustrates the sensitivity of fairness to the population model; in particular, none of the decision trees syntactically access
sex, yet several are unfair.
α
Figure 9 also shows that the modifications for affirmative action in DTα
16 and SVM 4 are sufficient to make the programs
fair with respect to all of our population models without making a substantial impact on the training set accuracy.
In summary, FairSquare is powerful enough to reason about group fairness for many non-trivial programs.
5.4

Comparison to other tools

We ran our benchmarks on the two other exact probabilistic inference tools that accept the same class of problems, and
report these results in the right bar graph in Figure 9. First, we compare to Sankaranarayanan et al.’s tool [66] (VC), which
13

2017/3/8

fairness ratio with sample
maximization

fairness ratio without sample
maximization

DT 16

SVM 4

4.0

4.0

3.5

DT16 nomax 5-step
DT16 nomax uniform
DT16 nomax none

3.0
2.5

(a)

2.0

3.5
2.5
1.5

1.0

1.0

0.5

0.5
0

5

10

15

20

25

30

4.0

0.0

0

5

10

15

20

25

30

4.0

3.5

DT16 max 5-step
DT16 max uniform
DT16 max none

3.0
2.5

(b)

2.0

3.5
2.5

1.0

1.0

0.5

0.5
5

10

15

20

25

(d)

2.0
1.5

0

SVM4 max 5-step
SVM4 max uniform
SVM4 max none

3.0

1.5

0.0

(c)

2.0

1.5

0.0

SVM4 nomax 5-step
SVM4 nomax uniform
SVM4 nomax none

3.0

30

0.0

0

5

10

15

20

25

30

Figure 10. Fairness ratio vs. rounds of sampling for DT16 and SVM4 (independent pop model) differing on ADFs and sample maximization. In (b) two
runs end at the exact value. Outside the visible range are: (a)(c) upper and lower bounds of uniform and none; (d) upper bounds of none.

is algorithmically similar to our tool: it finds bounds for probabilities on individual paths by approximating formulas with
bounding and inscribed hyperrectangles. Second, we compare to PSI [44], which symbolically computes representations
of the posterior distributions of variables. When closed form solutions to the CDFs exist, PSI is a great aid because
these solutions can be evaluated to obtain exact values of probabilities instantly; however, the graph illustrates that our
benchmarks are often too complex, resulting in integrals that do not have closed forms, or problems that PSI cannot solve
within the timeout period.
Tools were deemed to have failed on a benchmark when they timed out after a 900s period. Additionally, the bounds on
probability quantities from VC were often not tight enough for a proof. See Appendix D for the full quantitative table of
results.
The figure illustrates some qualitative properties of the applicability of the tools. In general, all can solve most of the
decision trees because the trees partition the decision space using a number of inequalities between a single variable and
a constant. However, the presence of inequalities involving multiple variables can result in two phenomena: (i) the lack
of closed form posterior CDFs, as reflected in the output of PSI, and (ii) angled boundaries in the decision space that are
hard to approximate with hyperrectangles. These inequalities occur in the SVMs, neural networks, and the Bayes Net 2
population model. Consequently, VC fails to produce good bounds in these cases; FairSquare performs better, but also fails
on the largest SVM.
5.5

Effect of Parameters

We now explore the effects of the approximate density functions (see Section 4.2) and of the sample maximization
optimizations. These results are captured in Figures 10 and 11.
There are three instances of ADFs in Figure 10 used to guide the sampling to high-probability regions: (i) none indicates
that no ADF is used, i.e., we used SYMVOL instead of ADF - SYMVOL; (ii) uniform indicates that each gauss(µ, σ 2 ) is
approximated by a uniform function spanning (µ − 3σ 2 , µ + 3σ 2 ) (similar to Figure 7(c)); and (iii) 5-step indicates
that each Gaussian is approximated by a step function of 5 equal-width regions spanning that same domain (similar to
Figure 7(a)). Another variable, max or nomax, denotes whether the sample maximization optimization is enabled (as
described in Section 5.1).
Each combination of these techniques is run on two of our benchmarks: DT16 and SVM4 under the independent
population model. Figure 10(a) and (c) show how convergence to the fairness ratio is improved by the choice of distribution
approximation when sample maximization is not employed: in particular, the runs using uniform and none are not even
visible, as the bounds never fall within [0.01, 4.0]. Plots (b) and (d) show that when sample maximization is employed, the
choice between the uniform and 5-step approximations is not as substantial on these benchmarks, although (i) the better
approximation gets better bounds faster, and (ii) using none results in substantially worse bounds.
14

2017/3/8

DT 16

average weighted
volume per sample

DT16

SVM 4
DT 16

nomax

0.15

SVM 4

max

SVM 4

nomax

0.15

(a)

(c)

0.1

0.1

0.05

0.05

0

0
5-step

average time per
round

max

uniform

none

5-step

2

uniform

none

2

(b)

1.5
1

1

0.5

0.5

0

(d)

1.5

0
5-step

uniform

none

5-step

uniform

none

Figure 11. Effect of optimizations on performance of FairSquare for DT16 and SVM4 (independent pop model). (a) and (c) show the average weighted
volume per sample (averaged across all computed probabilities). (b) and (d) show the average time (s) per round of sampling.

Figure 11 plot (a) and (c) show that employing ADFs and using sample maximization each increases the average
weighted volume per sample, which allows volume computation to be done with fewer samples. Plots (b) and (d) illustrate
the trade-off: the average time per round of sampling tends to be greater for more sophisticated optimizations.
We present these results for two particular problems and observe the same results across our suite. In summary, we have
found that ADFs and sample maximization are both necessary for adequate performance of FairSquare.

6.

Related Work

Probabilistic program analysis We refer the reader to Gordon et al. for a thorough survey [47]. A number of works tackled
analysis of probabilistic programs from an abstract interpretation perspective [29, 54, 57–59]. The comparison between
our solution through volume computation and abstract interpretation is perhaps analogous to SMT solving and software
model checking versus abstract interpretation. For example, techniques like Monniaux’s [57], sacrifice precision (through
joins, abstraction, etc.) of the analysis for the benefit of efficiency. Our approach, on the other hand, is aimed at eventually
producing a proof, or iteratively improving probability bounds while guaranteeing convergence.
A number of works have also addressed probabilistic analysis through symbolic execution [42, 45, 65, 66]. Filieri et
al. [42] and Geldenhuys et al. [45] attempt to find the probability a safety invariant is preserved. Both methods reduce
to a weighted model counting approach and are thus effectively restricted to variables over finite domains. Note that our
technique is more general than a model counting approach, as we can handle the discrete cases with a proper encoding of
the variables into a continuous domain.
A number of works rely on sampling to approximate the probability of a given program property. The Church [46]
programming language, for instance, employs the Metropolis–Hastings algorithm [26], a Markov Chain Monte Carlo
(MCMC) technique. Other techniques perform probabilistic inference by compiling programs or program paths to Bayesian
networks [51] and applying belief propagation [56] or sampling [65] on the network.
Volume computation The computation of weighted volume is known to be hard—even for a polytope, volume computation
is #P-hard [49]. Two general approaches exist: approximate and exact solutions. Note that in general, any approximate
technique at best can prove facts with high probability.
Our volume computation algorithm is inspired by (i) Li et al.’s [53] formula decomposition procedure, where quantifier
elimination is used to under-approximate an LRA constraint as a Boolean combination of monadic predicates; and
(ii) Sankaranarayanan et al.’s [66] technique for bounding the weighted volume of a polyhedron, which is the closest
volume computation work to ours. (The general technique of approximating complex regions with unions of orthogonal
polyhedra is well-studied in hybrid systems literature [21].)
A number of factors differentiate our work from [66], which we compared with experimentally in Section 5. First, our
approach is more general, in that it can operate on Boolean formulas over linear and polynomial inequalities, as opposed to
just conjunctions of linear inequalities. Second, our approach employs approximate distributions to guide the sampling of
hyperrectangles with large volume, which, as we have demonstrated experimentally, is a crucial feature of our approach.
Third, we provide theoretical convergence guarantees.
15

2017/3/8

LattE is a tool that performs exact integration of polynomial functions over polytopes [34]. Belle et al. [16, 18] compute
the volume of a linear real arithmetic (LRA) formula by, effectively, decomposing it into DNF—a set of polyhedra—and
using LattE to compute the volume of each polyhedron with respect to piece-wise polynomial densities. Our volume
computation algorithm is more general in that it (i) handles formulas over real closed fields, which subsumes LRA, and
(ii) handles arbitrary probability distributions.
Chistikov et al. [27, 28] present a framework for approximate counting with probabilistic guarantees in SMT theories,
which they specialize for bounded LRA. In contrast, our technique (i) handles unbounded formulas in LRA as well as real
closed fields, (ii) handles arbitrary distributions, and (iii) provides converging lower-bound guarantees. It is important to
note that there is a also a rich body of work investigating randomized polynomial algorithms for approximating the volume
of a polytope, beginning with Dyer et al.’s seminal work [38] (see Vempala [70] for a survey).
Algorithmic fairness Our work is inspired by recent concern in the fairness of modern decision-making programs [12, 71].
A number of recent works have explored algorithmic fairness [22, 32, 33, 37, 41, 48, 61, 72]. For instance, Zemel et al. [72]
and Feldman et al. [41] study fairness from a machine learning classification perspective, e.g., automatically learning fair
classifiers. Both works operate with a notion of fairness on a provided data set; in contrast, we prove fairness with respect
to a given probabilistic model of the population (a data set can be viewed as a special case).
Discrimination in black-box systems has been studied through the lens of statistical analysis [32, 33, 67]. Notably,
Datta et al. [32] created an automated tool that analyzes online advertising: it operates dynamically by surveying the ads
produced by Google. Our approach differs from the statistical analyses in that we require transparency of the decision
procedure, instead of data on the results of the decision procedure.
Future work We presented a new technique for verification of probabilistic programs and applied to quantifying bias
of decision-making programs. Our goals for future work include exploring the applicability of our technique to other
notions of fairness properties (e.g., [37]); further refinements to the use of ADFs, e.g. dynamically adapting them during
ADF - SYMVOL ; and improving the weighted volume computation algorithm for better scalability.

References
[1] Uci machine learning repository: Census income.
[2] Wolfram mathematica.

https://archive.ics.uci.edu/ml/datasets/Adult/.

https://www.wolfram.com/mathematica/,

[3] Matlab.

http://www.mathworks.com/products/matlab/,

[4] Redlog.

http://www.redlog.eu/.

[5] Weka.

.

.

http://www.cs.waikato.ac.nz/ml/weka/.

[6] Code of federal regulations.
https://www.gpo.gov/fdsys/pkg/CFR-2014-title29-vol4/xml/CFR-2014-title29-vol4part1607.xml, July 2014. (Accessed on 06/18/2016).
[7] I. Ajunwa, S. Friedler, C. E. Scheidegger, and S. Venkatasubramanian. Hiring by algorithm: predicting and preventing disparate
impact. Available at SSRN 2746078, 2016.
[8] J. Angwin, J. Larson, S. Mattu, and L. Kirchner. Machine bias: There’s software used across the country to predict future criminals. and it’s biased against blacks. https://www.propublica.org/article/machine-bias-risk-assessments-in-criminalsentencing, May 2016. (Accessed on 06/18/2016).
[9] E. Asarin, O. Bournez, T. Dang, and O. Maler. Approximate reachability analysis of piecewise-linear dynamical systems. In
International Workshop on Hybrid Systems: Computation and Control, pages 20–31. Springer, 2000.
[10] Y. Attali and J. Burstein. Automated essary scoring with e-rater v.2. Journal of Technology, Learning, and Assessment, 4(3), 2006.
[11] M. Barnett and K. R. M. Leino. Weakest-precondition of unstructured programs. In ACM SIGSOFT Software Engineering Notes,
volume 31, pages 82–87. ACM, 2005.
[12] S. Barocas and A. D. Selbst. Big data’s disparate impact. Available at SSRN 2477899, 2014.
[13] G. Barthe, P. R. D’Argenio, and T. Rezk. Secure information flow by self-composition. In CSFW, 2004.
[14] G. Barthe, M. Gaboardi, E. J. G. Arias, J. Hsu, C. Kunz, and P. Strub. Proving differential privacy in hoare logic. In IEEE
27th Computer Security Foundations Symposium, CSF 2014, Vienna, Austria, 19-22 July, 2014, pages 411–424, 2014. doi:
10.1109/CSF.2014.36. URL http://dx.doi.org/10.1109/CSF.2014.36.
[15] O. Bastani, Y. Ioannou, L. Lampropoulos, D. Vytiniotis, A. Nori, and A. Criminisi. Measuring neural net robustness with
constraints. CoRR, abs/1605.07262, 2016. URL http://arxiv.org/abs/1605.07262.
[16] V. Belle, A. Passerini, and G. V. den Broeck. Probabilistic inference in hybrid domains by weighted model integration. In
Proceedings of the Twenty-Fourth International Joint Conference on Artificial Intelligence, IJCAI 2015, Buenos Aires, Argentina,
July 25-31, 2015, pages 2770–2776, 2015. URL http://ijcai.org/Abstract/15/392.

16

2017/3/8

[17] V. Belle, G. Van den Broeck, and A. Passerini. Hashing-based approximate probabilistic inference in hybrid domains. In
Proceedings of the 31st Conference on Uncertainty in Artificial Intelligence (UAI), 2015.
[18] V. Belle, G. V. den Broeck, and A. Passerini. Component caching in hybrid domains with piecewise polynomial densities. In
Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, February 12-17, 2016, Phoenix, Arizona, USA., pages
3369–3375, 2016. URL http://www.aaai.org/ocs/index.php/AAAI/AAAI16/paper/view/12542.
[19] N. Berg. Predicting crime, lapd-style.

https://www.theguardian.com/cities/2014/jun/25/predicting-crime-lapd-los-

angeles-police-data-analysis-algorithm-minority-report,

June 2014. (Accessed on 06/18/2016).

[20] C. M. Bishop. Pattern recognition. Machine Learning, 128, 2006.
[21] O. Bournez, O. Maler, and A. Pnueli. Orthogonal Polyhedra: Representation and Computation, pages 46–60. Springer Berlin
Heidelberg, Berlin, Heidelberg, 1999. ISBN 978-3-540-48983-2. doi: 10.1007/3-540-48983-5_8. URL http://dx.doi.org/10.
1007/3-540-48983-5_8.
[22] T. Calders and S. Verwer. Three naive bayes approaches for discrimination-free classification. Data Mining and Knowledge
Discovery, 21(2):277–292, 2010.
[23] M. Carbin, S. Misailovic, and M. C. Rinard. Verifying quantitative reliability for programs that execute on unreliable hardware.
In Proceedings of the 2013 ACM SIGPLAN International Conference on Object Oriented Programming Systems Languages &
Applications, OOPSLA 2013, part of SPLASH 2013, Indianapolis, IN, USA, October 26-31, 2013, pages 33–52, 2013. doi:
10.1145/2509136.2509546. URL http://doi.acm.org/10.1145/2509136.2509546.
[24] S. Chaudhuri, S. Gulwani, R. Lublinerman, and S. Navidpour. Proving programs robust. In Proceedings of the 19th ACM SIGSOFT
Symposium and the 13th European Conference on Foundations of Software Engineering, ESEC/FSE ’11, pages 102–112, New
York, NY, USA, 2011. ACM. ISBN 978-1-4503-0443-6. doi: 10.1145/2025113.2025131. URL http://doi.acm.org/10.1145/
2025113.2025131.
[25] S. Chaudhuri, M. Clochard, and A. Solar-Lezama. Bridging boolean and quantitative synthesis using smoothed proof search. In
POPL, volume 49, pages 207–220. ACM, 2014.
[26] S. Chib and E. Greenberg. Understanding the metropolis-hastings algorithm. The american statistician, 49(4):327–335, 1995.
[27] D. Chistikov, R. Dimitrova, and R. Majumdar. Approximate counting in SMT and value estimation for probabilistic programs.
CoRR, abs/1411.0659, 2014. URL http://arxiv.org/abs/1411.0659.
[28] D. Chistikov, R. Dimitrova, and R. Majumdar. Approximate counting in SMT and value estimation for probabilistic programs. In
Tools and Algorithms for the Construction and Analysis of Systems - 21st International Conference, TACAS 2015, Held as Part of
the European Joint Conferences on Theory and Practice of Software, ETAPS 2015, London, UK, April 11-18, 2015. Proceedings,
pages 320–334, 2015. doi: 10.1007/978-3-662-46681-0_26. URL http://dx.doi.org/10.1007/978-3-662-46681-0_26.
[29] G. Claret, S. K. Rajamani, A. V. Nori, A. D. Gordon, and J. Borgström. Bayesian inference using data flow analysis. In Proceedings
of the 2013 9th Joint Meeting on Foundations of Software Engineering, pages 92–102. ACM, 2013.
[30] E. Clarke, D. Kroening, and F. Lerda. A tool for checking ansi-c programs. In International Conference on Tools and Algorithms
for the Construction and Analysis of Systems, pages 168–176. Springer, 2004.
[31] R. Cytron, J. Ferrante, B. K. Rosen, M. N. Wegman, and F. K. Zadeck. Efficiently computing static single assignment form and
the control dependence graph. ACM Transactions on Programming Languages and Systems (TOPLAS), 13(4):451–490, 1991.
[32] A. Datta, M. C. Tschantz, and A. Datta. Automated experiments on ad privacy settings. Proceedings on Privacy Enhancing
Technologies, 2015(1):92–112, 2015.
[33] A. Datta, S. Sen, and Y. Zick. Algorithmic transparency via quantitative input influence. In Proceedings of 37th IEEE Symposium
on Security and Privacy, 2016.
[34] J. De Loera, B. Dutra, M. Koeppe, S. Moreinis, G. Pinto, and J. Wu. Software for exact integration of polynomials over polyhedra.
ACM Communications in Computer Algebra, 45(3/4):169–172, 2012.
[35] L. De Moura and N. Bjørner. Z3: An efficient smt solver. In International conference on Tools and Algorithms for the Construction
and Analysis of Systems, pages 337–340. Springer, 2008.
[36] C. Dwork. Differential privacy. In Automata, languages and programming, pages 1–12. Springer, 2006.
[37] C. Dwork, M. Hardt, T. Pitassi, O. Reingold, and R. S. Zemel. Fairness through awareness. In Innovations in Theoretical Computer
Science 2012, Cambridge, MA, USA, January 8-10, 2012, pages 214–226, 2012. doi: 10.1145/2090236.2090255.
[38] M. Dyer, A. Frieze, and R. Kannan. A random polynomial-time algorithm for approximating the volume of convex bodies. Journal
of the ACM (JACM), 38(1):1–17, 1991.
[39] M. E. Dyer and A. M. Frieze. On the complexity of computing the volume of a polyhedron. SIAM Journal on Computing, 17(5):
967–974, 1988.
[40] V. Eubanks. The dangers of letting algorithms enforce policy.

http://www.slate.com/articles/technology/future_tense/

2015/04/the_dangers_of_letting_algorithms_enforce_policy.html,

17

April 2015. (Accessed on 06/18/2016).

2017/3/8

[41] M. Feldman, S. A. Friedler, J. Moeller, C. Scheidegger, and S. Venkatasubramanian. Certifying and removing disparate impact.
In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Sydney, NSW,
Australia, August 10-13, 2015, pages 259–268, 2015. doi: 10.1145/2783258.2783311. URL http://doi.acm.org/10.1145/
2783258.2783311.
[42] A. Filieri, C. S. Păsăreanu, and W. Visser. Reliability analysis in symbolic pathfinder. In Proceedings of the 2013 International
Conference on Software Engineering, pages 622–631. IEEE Press, 2013.
[43] S. A. Friedler, C. Scheidegger, and S. Venkatasubramanian. On the (im)possibility of fairness. CoRR, abs/1609.07236, 2016. URL
http://arxiv.org/abs/1609.07236.
[44] T. Gehr, S. Misailovic, and M. Vechev. Psi: Exact symbolic inference for probabilistic programs. In Computer aided verification.
Springer, 2016.
[45] J. Geldenhuys, M. B. Dwyer, and W. Visser. Probabilistic symbolic execution. In Proceedings of the 2012 International Symposium
on Software Testing and Analysis, pages 166–176. ACM, 2012.
[46] N. D. Goodman, V. K. Mansinghka, D. M. Roy, K. Bonawitz, and J. B. Tenenbaum. Church: a language for generative models. In
UAI 2008, Proceedings of the 24th Conference in Uncertainty in Artificial Intelligence, Helsinki, Finland, July 9-12, 2008, pages
220–229, 2008.
[47] A. D. Gordon, T. A. Henzinger, A. V. Nori, and S. K. Rajamani. Probabilistic programming. In Proceedings of the on Future of
Software Engineering, pages 167–181. ACM, 2014.
[48] M. Hardt, E. Price, and N. Srebro. Equality of opportunity in supervised learning. CoRR, abs/1610.02413, 2016. URL
http://arxiv.org/abs/1610.02413.
[49] L. Khachiyan. Complexity of polytope volume computation. Springer, 1993.
[50] N. Kobie. Who do you blame when an algorithm gets you fired?
accountable, January 2016. (Accessed on 06/18/2016).

http://www.wired.co.uk/article/make-algorithms-

[51] D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
[52] D. Kozen. Semantics of probabilistic programs. Journal of Computer and System Sciences, 22(3):328–350, 1981.
[53] Y. Li, T. H. Tan, and M. Chechik. Management of time requirements in component-based systems. In FM 2014: Formal Methods,
pages 399–415. Springer, 2014.
[54] P. Mardziel, S. Magill, M. Hicks, and M. Srivatsa. Dynamic enforcement of knowledge-based security policies. In Computer
Security Foundations Symposium (CSF), 2011 IEEE 24th, pages 114–128. IEEE, 2011.
[55] C. C. Miller. Can an algorithm hire better than a human? http://www.nytimes.com/2015/06/26/upshot/can-an-algorithmhire-better-than-a-human.html, June 2015. (Accessed on 06/18/2016).
[56] T. Minka, J. Winn, J. Guiver, and D. Knowles. Infer.net 2.5. Microsoft Research Cambridge, 2012.
[57] D. Monniaux. Abstract interpretation of probabilistic semantics. In Static Analysis, pages 322–339. Springer, 2000.
[58] D. Monniaux. An abstract monte-carlo method for the analysis of probabilistic programs. In ACM SIGPLAN Notices, volume 36,
pages 93–101. ACM, 2001.
[59] D. Monniaux. Backwards abstract interpretation of probabilistic programs. In Programming Languages and Systems, pages 367–
382. Springer, 2001.
[60] V. Nair and G. E. Hinton. Rectified linear units improve restricted boltzmann machines. In Proceedings of the 27th International
Conference on Machine Learning (ICML-10), pages 807–814, 2010.
[61] D. Pedreshi, S. Ruggieri, and F. Turini. Discrimination-aware data mining. In Proceedings of the 14th ACM SIGKDD international
conference on Knowledge discovery and data mining, pages 560–568. ACM, 2008.
[62] W. L. Perry. Predictive policing: The role of crime forecasting in law enforcement operations. Rand Corporation, 2013.
[63] J. Rawls. A theory of justice. Harvard university press, 2009.
[64] S. Ruggieri. Using t-closeness anonymity to control for non-discrimination. Transactions on Data Privacy, 7(2):99–129, 2014.
[65] A. Sampson, P. Panchekha, T. Mytkowicz, K. S. McKinley, D. Grossman, and L. Ceze. Expressing and verifying probabilistic
assertions. In ACM SIGPLAN Notices, volume 49, pages 112–122. ACM, 2014.
[66] S. Sankaranarayanan, A. Chakarov, and S. Gulwani. Static analysis for probabilistic programs: inferring whole program properties
from finitely many paths. In ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI ’13,
Seattle, WA, USA, June 16-19, 2013, pages 447–458, 2013. doi: 10.1145/2462156.2462179. URL http://doi.acm.org/10.
1145/2462156.2462179.
[67] L. Sweeney. Discrimination in online ad delivery. Queue, 11(3):10, 2013.
[68] A. Tutt. An fda for algorithms. Available at SSRN 2747994, 2016.

18

2017/3/8

s0 = s[x 7→ c]

s0 = s[x 7→ s(e)]
(ω, s, x ← e) → (ω, s0 , ε)

ASSIGN

(ω, s, S1 ) → (ω 0 , s0 , S10 )

(c : ω, s, x ∼ p) → (ω, s0 , ε)

PASSIGN

(ω, s, S1 S2 ) → (ω 0 , s0 , S10 S2 )

SEQ

s(b) = true
TCOND

SKIP

(ω, s, εS) → (ω, s, S)

(ω, s, if b then S1 else S2 ) → (ω, s, S1 )
s(b) = false
FCOND

(ω, s, if b then S1 else S2 ) → (ω, s, S2 )
Figure 12. Operational semantics

[69] J. Valentino-Devries, J. Singer-Vine, and A. Soltani. Websites vary prices, deals based on users’ information. http://www.wsj.
com/articles/SB10001424127887323777204578189391813881534, December 2012. (Accessed on 06/18/2016).
[70] S. Vempala. Geometric random walks: a survey. Combinatorial and computational geometry, 52(573-612):2, 2005.
[71] T. Zarsky. Understanding discrimination in the scored society. Washington Law Review, 89(4), 2014.
[72] R. S. Zemel, Y. Wu, K. Swersky, T. Pitassi, and C. Dwork. Learning fair representations. In Proceedings of the 30th
International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013, pages 325–333, 2013. URL
http://jmlr.org/proceedings/papers/v28/zemel13.html.

A.

Operational Semantics

Operational semantics We now define the operational semantics of our program model. Typically, the state s : V → R
of the program is defined as a valuation function from variables in V to values in R. In a probabilistic setting, however,
we need to maintain an additional state that dictates values drawn from probability distributions. Following standard
semantics of probabilistic programs [52], we assume a finite sequence of independent and identically distributed (iid)
random variables. The semantics of an execution is thus defined for a fixed sequence of values ω of these variables.
Informally, it is as if we performed all sampling before the program executes and stored the results in a sequence for use
whenever we encounter a probabilistic assignment. For simplicity of exposition, we assume all probabilistic assignments
sample from the same distribution; otherwise, multiple sequences can be used, one per distribution.
We can now define the semantics of P as shown in Figure 12. We use s(e) to denote the value of expression e under a
given state s, and use ε to denote the empty statement sequence. We use the substitution notation s[x 7→ c] to denote the
state s but with x mapped to c ∈ R. The interesting rule is PASSIGN, which, given a probabilistic assignment x ∼ p, for
some variable x and distribution p, picks and removes the first element c of the sequence c : ω to update the value of the
variable x.
We assume no variable is used before being assigned to. In a closed program, we assume all variables are initially
assigned to 0. Given a closed program P, we define JPKω as the final state s reachable from executing P using the
sequence ω, as defined by the relation → in Figure 12.

B.

Proofs

In this section, we prove correctness of the various pieces of our algorithm. First, we introduce preliminary theorems.
The following theorem states the soundness and completeness of hyperrectangular decomposition: models of ϕ
characterize all hyperrectangles in ϕ and no others.
Theorem 4 (Correctness of ). Let ϕ ∈ LT . Soundness: Let m |= ϕ . Then,VH m ⇒ ϕ is valid. Completeness: Let H
be a hyperrectangle such that H ⇒ ϕ. Then, the following is satisfiable: ϕ ∧ x∈Xϕ lx = Hl (x) ∧ ux = Hu (x)

The following theorem states the correctness of block : it removes all hyperrectangles that overlap with H m (soundness),
and it does not overconstrain Ψ by removing hyperrectangles that do not overlap with H m (completeness).

Theorem 5 (Correctness of block ). Given ϕ, let Ψ ⇒ ϕ , and let m1 , m2 |= Ψ. Soundness: If H m1 ∧ H m2 is satisfiable,
then m2 6|= Ψ ∧ block (H m1 ). Completeness: If H m1 ∧ H m2 is unsatisfiable, then m2 |= Ψ ∧ block (H m1 ).
The following theorem states the correctness of the ADF encoding:
Theorem 6 (Correctness of step φ ). Fix an ADF step(x). Soundness: For any model m |= step φ , the following is true:
R m(u )
Rb
m(δx ) = m(lxx) step(x) dx. Completeness: For any constants a, b, c ∈ R such that c = a step(x) dx, the following
formula is satisfiable: δx = c ∧ lx = a ∧ ux = b ∧ step φ (x).
19

2017/3/8

Soundness Suppose m |= ϕ . By definition of ϕ , the following formula is valid
Proof of Theorem 4

^
x∈Xϕ

m(lx ) 6 x 6 m(ux ) ⇒ ϕ

Therefore H m ⇒ ϕ, by definition of H m .
V
Completeness Let H be a hyperrectangle such that H ⇒ ϕ. By definition, H is of the form x∈Xϕ cx 6 x 6 c0x . It
immediately follows that the model where lx = cx and ux = c0x , for every x ∈ Xϕ , satisfies ϕ , since c0x > cx (satisfying
the first conjunct of ϕ ), and ∀Xϕ . H ⇒ ϕ (satisfying the second conjunct of ϕ ).
Proof of Theorem 5
Soundness Suppose H m1 ∧ H m2 is satisfiable. By definition of a hyperrectangle, this means that for all variables
x ∈ Xϕ , we have that the intervals [Hlm1 (x), Hum1 (x)] and [Hlm2 (x), Hum2 (x)] overlap, i.e., at least are equal on one
of the extremes. Therefore, m2 6|= Ψ ∧ block (H m1 ), since block does not admit any model m where, for all x ∈ Xϕ ,
[m(lx ), m(ux )] overlaps with [Hlm1 (x), Hum1 (x)].
Completeness Suppose H m1 ∧ H m2 is unsatisfiable. By definition of a hyperrectangle, there is at least one x ∈ Xϕ where
[Hlm1 (x), Hum1 (x)] and [Hlm2 (x), Hum2 (x)] do not overlap. Therefore, if m2 |= Ψ, then m2 |= Ψ ∧ block (H m1 ), since
block explicitly states that for at least one variable x, [m(lx ), m(ux )] should not overlap with [Hlm1 (x), Hum1 (x)].
Proof of Theorem 1
Pl R Q
At any point in the execution, vol = i=1 Hi pi (xi ) dxi , where l is the number of applications of HSAMPLE and Hi
W
is the hyperrectangle sampled at step i. By definition, Hi ⇒ ϕ. Since PDFs are positive functions, vol 6 VOL(ϕ, D).
Proof of Corollary 1
By definition of PDFs and integration,
Z Y

Z Y

pi (xi ) dxi =

Rn

Z
pi (xi ) dxi +

ϕ

Y

pi (xi ) dxi

¬ϕ

for any ϕ ⊆ Rn . From Theorem 1, it follows that at any point in the execution of SYMVOL(¬ϕ, D), we have 1 − vol >
VOL (ϕ, D).
Proof of Theorem 6
Soundness Suppose m |= step φ (x). Then,
m(δx ) =

n
X
i=1

ci · [ai , bi ] ∩ [m(lx ), m(ux )]

By definition of the area under a positive step function, we have
Z

m(ux )

m(δx ) =

step(x) dx
m(lx )

Completeness Completeness easily follows from correctness of the encoding of integrals over step functions as sums.
Proof of Theorem 2
P
The
vi and the approximated series
P algorithm constructs two series in parallel: the actual volume computation series
ai , where each vi and ai correspond to the actual and approximate volume of the i’th sampled hyperrectangle (note that
the latter is not explicitly maintained in the algorithm). Each series corresponds to a sequence of partial sums: Let
viΣ =

i
X

aΣ
i =

vj

j=1

i
X

aj

j=1

20

2017/3/8

It is maintained that
∀i. viΣ 6 EVolR∩ϕ =

Z

∀i. aΣ
i 6 AVol =

Z

Y

p(x) dXϕ

Y

step(x) dXϕ

R∩ϕ

R∩ϕ

Suppose, for the sake of obtaining a contradiction, that our sequence of samples to Rconstruct
{viΣ } and {aΣ
i } never, in
Q
0
0
the limit, samples some subregion R ( R, but H ⊆ R \ R is
a
hypercube
such
that
p(x)
dX
is
non-zero.
Then,
ϕ
R Q
RH Q
{viΣ } approaches some limit that is at most v Σ 6 EVolR∩ϕ − H p(x) dXϕ . Let δ = H step(x) dXϕ . The sequence
Σ
Σ
{aΣ
i } is monotonically increasing and bounded above by AVol: it follows that {ai } converges to some limit a 6 AVol.
Σ
Σ
By the definition of a limit, for all  > 0, there exists N such that for all n > N , a − an < . So at some point when we
have fixed a threshold τ < δ and have run out of samples in R0 with an > τ (guaranteed when aΣ − aΣ
n < τ ) we would
sample H ⊆ R \ R0 . This property ensures that in the limit, R0 → R, and viΣ → EVolR∩ϕ .
Proof of Theorem 3
More formally, ADF - SYMVOL using each Pi and its Ai infinitely often in the fair serialization translates to the existence
of a surjective function f : N → N such that (i) the nth successful call to HSAMPLE is in Pf (n) and (ii) for all n0 the
preimage f −1 (n0 ) is an infinite set.
S∞
Since the domain (Ri ) of each Ai is disjoint from the others, since i=1 Ri = Rn , and since each Pi converges to the
weighted volume over its domain, it follows that
∞ X
∞
X
i=1 j=1

VOL(Hij , D)

= VOL(ϕ, D)

where Hij is the jth hyperrectangle returned from Pi . (The nth hyperrectangle produced from the serialization corresponds
to Hij with i = f (n) and j = |f −1 (i) ∩ {n0 | n0 6 n}|.)
In fact, since all the terms are non-negative, the series above converges absolutely, and any rearrangement converges to
the same limit. A diagonalization argument permits a bijection b(i, j) = k such that this sum is equivalent to
∞
X
k=1

VOL(Hk , D)

= VOL(ϕ, D)

Thus, by composing rearrangements, any fair serialization also converges.

C.

Example Programs

As an example, below is the code for SVMα
4 with the independent population model. Note that the parameters for Gaussians
are mean and variance, and step functions are a list of tuples of the form (lower, upper, value). A decision-making program
F is free to use variables defined in the population model popModel.
def popModel(): //independent
age ~ gaussian(38.5816, 186.0614)
sex ~ step([(0,1,0.3307), (1,2,0.6693)])
capital_gain ~ gaussian(1077.6488, 54542539.1784)
capital_loss ~ gaussian(87.3038, 162376.9378)
sensitiveAttribute(sex < 1)
qualified(True)
def F(): //SVM4A
N_age = (age - 17.0) / 62.0
N_capital_gain = capital_gain / 22040.0
N_capital_loss = capital_loss / 1258.0
t = 0.0006 * N_age
+ -5.7363 * N_capital_gain
+ -0.0002 * N_capital_loss
+ 1.0003
if sex > 1:
t = t + -0.0003

21

2017/3/8

if sex < 1: //affirmative action
t = t - 0.5
fairnessTarget(t < 0)

Below is the Bayes Net 2 population model (Bayes Net 1 is the same without last if statement):
def popModel(): //Bayes Net 2
sex ~ step([(0,1,0.3307), (1,2,0.6693)])
if sex < 1:
capital_gain ~ gaussian(568.4105, 24248365.5428)
if capital_gain < 7298.0000:
age ~ gaussian(38.4208, 184.9151)
education_num ~ gaussian(10.0827, 6.5096)
capital_loss ~ gaussian(86.5949, 157731.9553)
else:
age ~ gaussian(38.8125, 193.4918)
education_num ~ gaussian(10.1041, 6.1522)
capital_loss ~ gaussian(117.8083, 252612.0300)
else:
capital_gain ~ gaussian(1329.3700, 69327473.1006)
if capital_gain < 5178.0000:
age ~ gaussian(38.6361, 187.2435)
education_num ~ gaussian(10.0817, 6.4841)
capital_loss ~ gaussian(87.0152, 161032.4157)
else:
age ~ gaussian(38.2668, 187.2747)
education_num ~ gaussian(10.0974, 7.1793)
capital_loss ~ gaussian(101.7672, 189798.1926)
if (education_num > age):
age = education_num
sensitiveAttribute(sex < 1)
qualified(True)

D.

Quantitative results

The full quantitative results of our experiments (Section 5) are shown in table 1; note that the upper-left quadrant of table 1
is simply the left table from figure 9. The additional rows are the instances of qualified group fairness where qualification
is true if and only if age > 18 (as opposed to the tautologically true case for qualification). The comparison to other tools
was run only on the trivial qualification, since the runs on FairSquare proved to be qualitatively very similar.
The second set of columns detail the results of running Sankaranarayanan et al.’s tool [66], VC. For the Res column,
a green check mark indicates that it proved the program fair, and a red strike indicates that it proved the program unfair.
Sometimes it would terminate with bounds that were not conclusive: the upper and lower bounds for the fairness ratio are
shown in Res. T denotes the amount of time (s) the tool ran for; TO denotes that it timed out after 900s without producing
a result.
The third set of columns detail the results of running PSI on the benchmarks. T denotes the amount of time (s) the tool
ran for; again, TO denotes that it timed out after 900s without producing a result. The Res column depicts a check mark if
PSI terminated with a closed form CDF , and an integral sign if it returned a function containing unevaluated integrals.

22

2017/3/8

FairSquare

VC

[66]

Decision
program

Acc
Res

#

Vol

QE

Res

#

Vol

QE

Res

#

Vol

QE

Res

T

Res

T

DT4

0.79

3

10

1.3

0.5

7

12

2.2

0.9

7

18

6.6

2.2

3

4.8

7

3.7

DT14

0.71

3

20

4.2

1.4

3

38

52.3

11.4

3

73 130.9 33.6

0.76
1.34

12.3

0.74
1.05

24.2

DT 16

0.79

3

21

7.7

2.0

7

22

15.3

6.3

7

22

38.2

14.3

3

18.7

7

22.3

DT16

0.76

3

18

5.1

3.0

3

34

32.0

8.2

3

40

91.0

19.4

3

23.2

3

36.0

DT44

0.82

3

55 63.5 9.8

7

113 178.9 94.3

7

406 484.0 222.4

3

40.3

SVM 3
SVM 4
α
SVM 4
SVM 5

0.79
0.79
0.78
0.79

3
3
3
3

10
10
10
10

7
7
3

10
14
14
-

10.8
33.7
33.4
-

6.2
20.1
63.2

TO
TO
TO
TO

295.8

TO

-

0.56
1.00
0.00
153.8

-

TO
TO
TO

SVM 6

0.79

0.02
35.3

-

-

TO

-

TO

-

TO

NN 2,1

0.65

3

78 21.6 0.8

3

466 456.1

3.4

3

154 132.9

NN 2,2

0.67

3

62 27.8 2.0

3

238 236.5

7.2

3

174 233.5 18.2

NN 3,2
DT4 Q
DT14 Q
DT16 Q
α
DT16 Q
DT44 Q
SVM 3 Q
SVM 4 Q

0.74

0.03
674.7

α

Independent

Bayes Net 1

2.6
2.7
3.0
8.5

0.6
0.8
0.8
1.3

7
7
3
7

10
18
22
10

3.7
13.3
15.7
12.2

634 TO

2.4

0.09
3.03

434

TO

0.00
5.24

3
3
3
3
3
3
3

442
12
24
27
16
50
22
18

2.0
5.7
10.1
5.3
67.5
3.1
4.4

10.0
1.1
2.4
3.6
3.7
13.1
3.0
3.9

SVM 4 Q
SVM 5 Q

3
3

14
18

4.2
6.4

1.9
5.3

SVM 6 Q

0.24
31.5

NN 2,1 Q
NN 2,2 Q

3
3

50 11.5 1.2
50 27.4 3.1

3
3

NN 3,2 Q

0.25
47.4

394 TO 14.3

0.07
1010

α

TO

686 TO

3.6

7
3
7
3
7
7
7
0.82
1.12

7
0.04
3.04

34
25
47
38
45
271
7.9
31

Bayes Net 2

1.7
3.1
3.2
6.3

TO q

12.8 TOq

7.2

23
80
128
340
24
24

10.4
57.1
179.7
337.9
18.3
34.8

705 TO
39 83.2

5.9
3 48
54.0 TOq -

94.9
-

49.4

440

133.4 TOq

-

TO

TO

5.6
57.1
48.9
46.2
347.5
27
42.9

TO

55.9 TOq
1.8
7
13.8 3
15.3 7
9.3
3
285.9 TOq
4.8
7
19.6 7

Ind

107 57.5 3.4
584 799.6 13.6
93

TO

0.71
1.09

3

53.4 TOq

-

0.00
691.3
0.00
1365

TO

-

PSI

BN1

BN2

56.3

Res
0.07
3.82
0.18
7.74
0.07
5.06
0.13
6.00
0.10
6.02

-

97.1

0.00
410.1
0.00
162.7

65.44

0.00
1012
0.00
67.1
0.01
29.2

TO

-

TO

-

27.4

42.3

Ind

[44]

BN1

T

Res

T

Res

T

16.9

4

0.07

4

0.22

81.2

6.7

BN2
Res
R

T
1.8

4

0.16

4

78.45 4

0.2

4 136.4

-

TO

4

0.6

-

TO

-

TO

TO

0.71
1.9
1.9
5.0

R
R
R
R

4.2
8.2
8.3
18.2

57.8

R

30.1

5.6

R

22.4

712.2

-

TO

TO

-

TO

99.8

214.6 4 539.3
TO
TO
TO
TO

4
4
4
4

0.06
0.09
0.09
0.14

R
R
R
R

R
248.0 4 0.24
R
R
80.7
1.1
R
R
107.1
541.8
R
TO
663.7 -

-

TO

TO

4.2
13.8
32.6
44.1
TO

7.4
148.3
TO

710 TO 11.8
431 558.5 21.5
-

-

TO

Table 1. Benchmarks

23

2017/3/8

