Quantifying Program Bias
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. %When we automate decision-making tasks, %we risk instilling bias in programs. We propose a novel probabilistic program analysis technique and apply it to quantifying bias in decision-making programs. Specifically, we \rone present a sound and complete automated verification technique for proving quantitative properties of probabilistic programs; \rtwo show that certain notions of bias, recently proposed in the fairness literature, can be phrased as quantitative correctness properties; and \rthree present \fs, the first verification tool for quantifying program bias, and evaluate it on a range of decision-making programs. % generated by a range of machine-learning % algorithms from real-world data. % Algorithms have become powerful % arbitrators of an ever-expanding range of sensitive % decisions with significant societal % impact---making hiring decisions, % deciding credit scores, deciding welfare eligibility, % predicting crime activity, directing policing, % influencing public policy, % amongst many others. % And the indicators are all clear: Algorithmic decision-making % is unlikely to stop in the foreseeable future; % in fact, we are merely witnessing its very beginning. % % With the range and sensitivity of algorithmic decisions % expanding at a break-neck speed, % it is imperative that we aggressively investigate % algorithmic fairness % and bias. % In this paper, we ask the following question: % Does a decision-making program discriminate % against a protected group? % Specifically, % we \rone formalize algorithmic fairness as a % probabilistic--program--verification problem, % and \rtwo present a novel sound and complete % automated verification technique, based on % volume computation, for % proving (un)fairness of decision-making programs. % % We implemented our techniques in a new fairness % verification tool called \fs, and applied % it to proving fairness of decision-making % programs generated by a range of machine learning % algorithms from real-world data. % Our results highlight the power and promise % of our approach at attacking this new and % increasingly important problem.
alg.tex
\section{Fairness-Verification Algorithm}\label{sec:alg}

\todo{Turn this into a case study.
Begin by introducing the notion of group fairness
similar to the current intro section.}

Now that we have fully introduced our volume-computation
technique, we are ready to present the fairness-verification algorithm.
Recall that in an $\epsilon$-fairness problem, $(\pre, \dec, \post)$,
we would like to prove or disprove the probabilistic postcondition $(\post)$:
$$\frac{\pr{\varphi_f \mid \varphi_s}}{
        \pr{\varphi_f \mid \neg \varphi_s}} > 1 - \epsilon$$
To prove this statement, we need to compute
a value for the numerator and the denominator. However,
recall that $\volalg$ maintains a \emph{lower bound} on the
probability. Therefore, we have to be careful in how we use its results.
Specifically, to prove fairness, we need a lower bound on
the numerator and an upper bound on the denominator---these bounds allow
us to compute a lower bound on the ratio.
Conversely, to prove unfairness, we require
an upper bound on the numerator
and a lower bound on the denominator.

First, let us rewrite $\post$ to get rid of conditional probabilities,
using the axiom $\pr{A \mid B} = \pr{A \land B} / \pr{B}$:
$$\frac{\pr{\varphi_f \land \varphi_s} \cdot \pr{\neg \varphi_s}}
{\pr{\varphi_f \land \neg \varphi_s} \cdot \pr{\varphi_s}}
> 1-\epsilon
$$
As in the abstract verification algorithm from Figure~\ref{alg:verify},
we have the probabilistic verification condition, $\pair{\varphi_\prog, \pdfs}$,
and the set of variables $V_d = V_d^\emph{pre} \cup V_d^\emph{dec}$.
To find a lower bound on the probability of some property $\varphi$,
denoted by $\prl{\varphi}$, we can invoke $\volalg(\exists V_d \ldotp \semprog \land \varphi, \pdfs)$.
Conversely, following Corollary~\ref{col:upper}, to find an upper bound $\pru{\varphi}$,
we can compute $1 - \volalg(\neg \exists V_d \ldotp \semprog \land \varphi, \pdfs).$

Our fairness algorithm is shown in Figure~\ref{alg:fairsquare}.
We assume that the function \textsc{updateProbabilities}
runs \volalg for a limited amount of iterations to
improve the upper and lower bounds on the various
properties, $\pr{\varphi_f \land \varphi_s}$, $\pr{\varphi_s}$, etc.
Thus, the algorithm simply loops, and in every iteration checks
whether we can prove fairness or unfairness.
Specifically, it first computes a lower bound
on the ratio
$\frac{\pr{\varphi_f \mid \varphi_s}}{
        \pr{\varphi_f \mid \neg \varphi_s}}$;
if the lower bound is greater than $1-\epsilon$,
then the problem is definitely fair.
The algorithm then proceeds to compute an upper
bound on the ratio, and if the upper bound is $\leq 1-\epsilon$,
we know that the problem is definitely unfair.

\paragraph{Soundness and completeness}
The following theorem states
soundness and (\emph{almost}) completeness
of \textsc{fairVerify}.
Completeness of the algorithm
is contingent on two factors.
First, \volalg should in the limit
converge to exact solutions.
If the problem is fair, then, in finitely many steps,
 we will be able to prove that the ratio is $> 1-\epsilon$.
However, in the rare chance that
the ratio $\frac{\pr{\varphi_f \mid \varphi_s}}{
        \pr{\varphi_f \mid \neg \varphi_s}}$
        is \emph{exactly}  $1 - \epsilon$,
we might not be able to prove that the algorithm
is unfair, as our upper bound on the ratio will come closer
and closer to $1 - \epsilon$, but potentially
never reach it.


\begin{theorem}[Correctness of \textsc{fairVerify}]\label{thm:fair}
Fix a fairness verification problem $(\pre, \dec, \post)$.
%
%\begin{itemize}
%
%\item
\textbf{Soundness}:
If \textsc{fairVerify} returns fair or unfair,
the result is correct.
%
%\item
\textbf{Almost completeness}:
If $\frac{\pr{\varphi_f \mid \varphi_s}}{
        \pr{\varphi_f \mid \neg \varphi_s}} \neq 1 - \epsilon$,
then \textsc{fairVerify} eventually
terminates.
%\end{itemize}
\end{theorem}

\begin{figure}[!t]
  \small
  \begin{spacing}{1}
    \begin{algorithmic}
        \Function{fairVerify}{$\pre, \dec, \post$}
          % \State Let $\varphi_\prog$ and $\pdfs$ be as in Figure~\ref{alg:verify}
          % \State Fork the following \volalg processes
          % \State $\prl{\varphi_f\land\varphi_s} \gets \volalg( \exists V_d \ldotp \varphi_\prog \land \varphi_s \land \varphi_f, \pdfs) $
          % \State $\pru{\varphi_f\land\varphi_s} \gets 1 - \volalg(\neg \exists V_d \ldotp \varphi_\prog \land \varphi_s \land \varphi_f, \pdfs) $
          % \State
          % \State $\prl{\varphi_f\land \neg \varphi_s} \gets \volalg( \exists V_d \ldotp \varphi_\prog \land \neg \varphi_s \land \varphi_f, \pdfs) $
          % \State $\pru{\varphi_f\land \neg\varphi_s} \gets 1 - \volalg(\neg \exists V_d \ldotp \varphi_\prog \land \neg \varphi_s \land \varphi_f, \pdfs) $
          % \State
          % \State $\prl{\varphi_s} \gets \volalg( \exists V_d \ldotp \varphi_\prog \land \varphi_s, \pdfs) $
          % \State $\pru{\varphi_s} \gets 1 - \volalg(\neg \exists V_d \ldotp \varphi_\prog \land \neg \varphi_s \land \varphi_f, \pdfs) $
          \While {$\true$}

			~~~~~~~\Call{updateProbabilities}{\hspace{0mm}}

            {\color{darkgray}{~~~~~~~// compute lower bound on ratio}}
            \begin{align*}
              \textbf{if}~~
            \frac{\prl{\varphi_f \land \varphi_s} \cdot \prl{\neg \varphi_s}}
            {\pru{\varphi_f \land \neg \varphi_s} \cdot \pru{\varphi_s}}
            > 1-\epsilon & &
          \end{align*}
               \hspace{.7in}\Return fair


               {\color{darkgray}{~~~~~~~// compute upper bound on ratio}}
              \begin{align*}
                \textbf{if}~~
              \frac{\pru{\varphi_f \land \varphi_s} \cdot \pru{\neg \varphi_s}}
              {\prl{\varphi_f \land \neg \varphi_s} \cdot \prl{\varphi_s}}
              \leq  1-\epsilon & &
            \end{align*}
                 \hspace{.7in}\Return unfair
            \EndWhile
        \EndFunction
    \end{algorithmic}
    \caption{Fairness verification algorithm}\label{alg:fairsquare}
  \end{spacing}
\end{figure}

bigtable.tex
\begin{table*}[]
\newcommand{\qentry}[9]{#1 & #2 & #3 & #4 & #5 & #6 & #7 & #8 & #9}
\centering
\begin{tabular}{|c|c|c|c|c||cc||cc|cc|cc|}
\cline{2-13}
\multicolumn{1}{c|}{} & Dist Approx & Sample Max & Round & Time (s) & \multicolumn{2}{c||}{Fairness} & \multicolumn{2}{c|}{$\pr{\varphi_s}$} & \multicolumn{2}{c|}{$\pr{\varphi_f \land \varphi_s}$} & \multicolumn{2}{c|}{$\pr{\varphi_f \land \lnot \varphi_s}$} \\ \hline
\multirow{24}{*}{DT4} & \multirow{4}{*}{none} & \multirow{4}{*}{no} 
                     & 1 & \qentry{0.56}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{1.00} \\
                 & & & 2 & \qentry{0.67}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{0.99} \\ 
                 & & & 3 & \qentry{0.78}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-12}$}{0.99} \\ 
                 & & & 4 & \qentry{0.90}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-12}$}{0.99} \\ \cline{2-13}
                     & \multirow{4}{*}{uniform} & \multirow{4}{*}{no} 
                     & 1 & \qentry{0.41}{0.01}{40.32}{0.33}{0.33}{0.00}{0.33}{0.02}{0.67} \\
                 & & & 2 & \qentry{0.57}{0.06}{13.23}{0.33}{0.33}{0.01}{0.33}{0.05}{0.34} \\ 
                 & & & 3 & \qentry{0.72}{0.08}{12.69}{0.33}{0.33}{0.01}{0.32}{0.05}{0.34} \\ 
                 & & & 4 & \qentry{0.88}{0.08}{7.60}{0.33}{0.33}{0.01}{0.29}{0.08}{0.34} \\ \cline{2-13}
                     & \multirow{4}{*}{5 steps} & \multirow{4}{*}{no} 
                     & 1 & \qentry{0.70}{0.21}{7.50}{0.33}{0.33}{0.05}{0.35}{0.09}{0.50} \\
                 & & & 2 & \qentry{1.01}{0.27}{3.98}{0.33}{0.33}{0.06}{0.23}{0.12}{0.41} \\ 
                 & & & 3 & \qentry{1.37}{0.43}{2.56}{0.33}{0.33}{0.06}{0.16}{0.12}{0.27} \\ 
                 & & & 4 & \qentry{1.80}{0.45}{2.20}{0.33}{0.33}{0.06}{0.14}{0.13}{0.26} \\ \cline{2-13}
                     & \multirow{4}{*}{none} & \multirow{4}{*}{yes} 
                     & 1 & \qentry{0.73}{0.00}{11.20}{0.33}{1.00}{0.08}{1.00}{0.18}{0.58} \\
                 & & & 2 & \qentry{0.90}{0.62}{11.20}{0.33}{0.33}{0.08}{1.00}{0.18}{0.29} \\ 
                 & & & 3 & \qentry{1.04}{0.63}{8.86}{0.33}{0.33}{0.09}{0.79}{0.18}{0.29} \\ 
                 & & & 4 & \qentry{1.20}{0.63}{8.86}{0.33}{0.33}{0.09}{0.79}{0.18}{0.29} \\ \cline{2-13}
                     & \multirow{4}{*}{uniform} & \multirow{4}{*}{yes} 
                     & 1 & \qentry{0.70}{0.27}{3.71}{0.33}{0.33}{0.09}{0.33}{0.18}{0.67} \\
                 & & & 2 & \qentry{0.86}{0.54}{2.00}{0.33}{0.33}{0.09}{0.18}{0.18}{0.33} \\ 
                 & & & 3 & \qentry{1.04}{1.00}{1.00}{0.33}{0.33}{0.09}{0.09}{0.18}{0.18} \\ 
                 & & & 4 & \qentry{1.14}{1.00}{1.00}{0.33}{0.33}{0.09}{0.09}{0.18}{0.18} \\ \cline{2-13}
                     & \multirow{4}{*}{5 steps} & \multirow{4}{*}{yes} 
                     & 1 & \qentry{0.98}{0.26}{3.70}{0.33}{0.33}{0.08}{0.33}{0.18}{0.66} \\
                 & & & 2 & \qentry{1.17}{0.50}{2.00}{0.33}{0.33}{0.09}{0.18}{0.18}{0.36} \\ 
                 & & & 3 & \qentry{1.36}{1.00}{1.00}{0.33}{0.33}{0.09}{0.09}{0.18}{0.18} \\ 
                 & & & 4 & \qentry{1.47}{1.00}{1.00}{0.33}{0.33}{0.09}{0.09}{0.18}{0.18} \\ \hhline{|=|=|=|=|=#==#==|==|==|}
\multirow{24}{*}{SVM3} & \multirow{4}{*}{none} & \multirow{4}{*}{no} 
                      & 1 & \qentry{0.57}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{1.00} \\
                  & & & 2 & \qentry{0.69}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{1.00} \\ 
                  & & & 3 & \qentry{0.80}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{1.00} \\ 
                  & & & 4 & \qentry{0.91}{0.00}{$10^{11}$}{0.33}{1.00}{$10^{-12}$}{1.00}{$10^{-11}$}{1.00} \\ \cline{2-13}
                      & \multirow{4}{*}{uniform} & \multirow{4}{*}{no} 
                      & 1 & \qentry{0.44}{0.02}{129.07}{0.25}{0.33}{0.01}{0.33}{0.01}{0.99} \\
                  & & & 2 & \qentry{0.65}{0.07}{19.46}{0.33}{0.33}{0.03}{0.33}{0.03}{0.98} \\ 
                  & & & 3 & \qentry{0.81}{0.11}{11.84}{0.33}{0.33}{0.04}{0.33}{0.06}{0.71} \\ 
                  & & & 4 & \qentry{0.99}{0.11}{11.01}{0.33}{0.33}{0.04}{0.32}{0.06}{0.70} \\ \cline{2-13}
                      & \multirow{4}{*}{5 steps} & \multirow{4}{*}{no} 
                      & 1 & \qentry{0.74}{0.09}{3.01}{0.25}{0.33}{0.03}{0.33}{0.33}{0.75} \\
                  & & & 2 & \qentry{1.31}{0.18}{1.15}{0.33}{0.33}{0.04}{0.19}{0.33}{0.42} \\ 
                  & & & 3 & \qentry{1.84}{0.27}{1.09}{0.33}{0.33}{0.06}{0.18}{0.33}{0.41} \\ 
                  & & & 4 & \qentry{2.45}{0.31}{1.08}{0.33}{0.33}{0.06}{0.18}{0.33}{0.39} \\ \cline{2-13}
                      & \multirow{4}{*}{none} & \multirow{4}{*}{yes} 
                      & 1 & \qentry{0.77}{0.00}{10.39}{0.33}{1.00}{0.10}{1.00}{0.19}{0.79} \\
                  & & & 2 & \qentry{1.10}{0.25}{8.21}{0.33}{0.33}{0.10}{0.79}{0.20}{0.79} \\ 
                  & & & 3 & \qentry{1.43}{0.25}{8.20}{0.33}{0.33}{0.10}{0.79}{0.20}{0.79} \\ 
                  & & & 4 & \qentry{1.83}{0.25}{8.19}{0.33}{0.33}{0.10}{0.79}{0.20}{0.79} \\ \cline{2-13}
                      & \multirow{4}{*}{uniform} & \multirow{4}{*}{yes} 
                      & 1 & \qentry{0.65}{0.34}{2.00}{0.33}{0.33}{0.17}{0.33}{0.33}{0.98} \\
                  & & & 2 & \qentry{1.03}{0.40}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.83} \\ 
                  & & & 3 & \qentry{1.47}{0.49}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.68} \\ 
                  & & & 4 & \qentry{1.99}{0.97}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.34} \\ \cline{2-13}
                      & \multirow{4}{*}{5 steps} & \multirow{4}{*}{yes} 
                      & 1 & \qentry{0.94}{0.25}{2.00}{0.33}{0.33}{0.09}{0.33}{0.33}{0.75} \\
                  & & & 2 & \qentry{1.51}{0.89}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.38} \\ 
                  & & & 3 & \qentry{2.38}{0.98}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.34} \\ 
                  & & & 4 & \qentry{3.44}{0.99}{1.00}{0.33}{0.33}{0.17}{0.17}{0.33}{0.34} \\ \hline
\end{tabular}
\caption{A table illustrating how each of our various optimizations affect
the performance of \fs.  The population model used throughout is
the set of independent variables.  Each instance shows the intermediate calculations
for each of the first 4 rounds of volume computation.
The four leftmost columns show the lower and upper bounds for their
corresponding values.}
\label{tbl:opt}
\end{table*}



code.tex
\section{Example Programs} \label{app:code}
As an example, below is the code for
\svma{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 \texttt{F} is free to use
variables defined in the population model \texttt{popModel}.

\begin{lstlisting}
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
    if sex < 1: //affirmative action
        t = t - 0.5
    fairnessTarget(t < 0)
\end{lstlisting}
Below is the Bayes Net 2 population model
(Bayes Net 1 is the same without last \texttt{if} statement):
\begin{lstlisting}
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)
\end{lstlisting}

conclusion.tex
\section{Conclusion}\label{sec:conclusion}

We addressed the verification problem of group fairness
and implemented our tool, \fs, 
that utilizes a novel weighted volume computation algorithm 
to prove the fairness or unfairness of programs.
%
Our goals for future work include
exploring the applicability of our technique
to \emph{individual fairness}~\cite{dwork12};
further refinements to the use of \abrs{ADF},
e.g. dynamically adapting them during \volopt;
and improving the weighted volume computation algorithm
for better scalability in the number of variables.


discussion.tex
\section{Discussion}
\label{sec:discussion}

In this section we discuss rationale behind the important design decisions that were made in this paper.

\todo{This paragraph is no longer very interesting}
We start by explaining why our new definition of group fairness uses a population model
instead of a concrete dataset.
Existing definitions of fairness define programs as fair or unfair with respect to a given concrete dataset, e.g.,~\cite{feldman15,dwork12}.
While using a concrete dataset simplifies the verification problem,
it also raises questions of whether
the dataset is representative for the population for which we are trying to prove fairness.
%Moreover, sensitive datasets are unlikely to be publicly available.
Our technique moves away from concrete datasets and replaces them with a  probabilistic population model.
Population models subsume concrete datasets and do not require making datasets publicly available.
Building a representative population model can be challenging; however, we envision a future in which
governmental agencies will publish representative population models and use them to assess
 fairness of algorithmic processes. The European Union (EU), for instance, has already begun regulating algorithmic decision-making~\cite{Goodman16}.

We presented an algorithm for exact volume computation (probabilistic inference) instead of using
existing techniques for approximate volume computation
or probabilistic inference, e.g. \emph{Metropolis--Hastings}~\cite{chib1995understanding} or
convex volume approximation~\cite{vempala2005geometric}.
Approximation techniques trade soundness for performance,
but soundness is fundamental for checking algorithm fairness:
%We can't conclude that an algorithm is more or less fair.
Our technique allows us to generate a checkable
\emph{certificate of fairness}, which is the set of sampled hyperrectangles.
Moreover, since our algorithm computes several individual quantities and then combines them,
inexact methods might end up accumulating error and producing very inaccurate results.

\todo{Need to rethink this justification}
While the technique we proposed is fairly general and captures many interesting machine-learning
classifiers, our goal is to build a fairness verifier and not a general-purpose
probabilistic program verifier.
For this reason, our language only contains constructs
that are used in classes of programs
where we might be interested in reasoning about fairness.
For instance, a majority of classifiers
synthesized by machine-learning algorithms are loop-free
programs that operate over real-valued inputs,
which our program model is able to capture.
% For example, our tool only allows distributions with constant parameters---e.g., mean and standard deviation.

%!TEX root = paper.tex

example.tex
\section{Overview and Illustration}\label{sec:example}

%\subsection{Problem Setting}
Our problem setting is as follows:
First, we are given a \emph{decision-making program} $\prog_\emph{dec}$.
Second, we have a \emph{probabilistic precondition} defining a probability
distribution over inputs of $\prog_\emph{dec}$.
We define the probability distribution operationally as a probabilistic
program $\prog_\emph{pre}$, which we call the \emph{population model}.
Intuitively, the population model provides a probabilistic picture
of the population from which the inputs of  $\prog_\emph{dec}$ are drawn.
% The probabilistic language we define
% is similar
% in expressiveness to recent works on probabilistic programming~\cite{nori2014r2,sankaranarayanan13,Gehr16,sampson2014expressing}.
Third, we are given a quantitative postcondition $\post$ that correlates the probabilities of various program outcomes.
%While we provide a general language for postconditions,
%the design of our language is inspired by
%recent bias and fairness definitions from the literature.
This postcondition can encode various properties
relating to program bias;
%
intuitively, our goal is to prove the following triple:
$${\cond{\vec{v} \sim \prog_\emph{pre}}}\quad
r \gets \prog_\emph{dec}(\vec{v})\quad
{\bcond{\post
}}$$

Let us consider various possible instantiations of $\post$.
Feldman et al.~\cite{feldman15} introduced the following
definition, inspired by \emph{Equality of Employment
Opportunity Commission's}~\cite{EEOC} recommendation
in the US:
\[
   \frac{\pr{r = \true  \mid \emph{min}(\vec{v}) = \true}}
    {\pr{r = \true  \mid \emph{min}(\vec{v}) = \false}}
>
1 - \epsilon
\]
Assuming $\prog_\emph{dec}$ returns a Boolean value---indicating
whether an applicant $\vec{v}$ is hired---%
this \emph{group fairness} property states that the selection rate
from a minority group, $\emph{min}(\vec{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
\emph{two sets of program traces}: one set where the input $\emph{min}(\vec{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 (\emph{statistical parity}).
Further, we could additionally condition on \emph{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., ~\cite{feldman15,zemel13,datta2016algorithmic}.

While the above definition is concerned with fairness at the level
of subsets of the domain of the decision-making program,
\emph{individual fairness}~\cite{dwork12} is concerned with similar outcomes
for similar elements of the domain.
In our hiring example, one potential formulation is as follows:
\[
\pr{r_1 = r_2  \mid \vec{v}_1 \sim \vec{v}_2 } > 1 - \epsilon
\]
In other words, we want to ensure that for any two
individuals, if they are \emph{similar} ($\sim$),
then we want them to receive similar outcomes ($r_1 = r_2$) with a high probability.
This is a \emph{hyperproperty}---as it considers two copies
of $\prog_\emph{dec}$---and can be encoded
through \emph{self-composition}~\cite{barthe04}.
This property is close in nature to differential privacy~\cite{dwork2006differential}
and robustness~\cite{Chaudhuri11,bastani16}.

Of course, various definitions of fairness have
their merits and their shortcomings,
and there is an ongoing discussion on this subject~\cite{friedler16,dwork12,hardt16,feldman15,ajunwa2016hiring}.
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.
% It is important to note that, while
% recent focus on program fairness in the literature
% and the media has been on automated bureaucratic
% processes that may cause social harm,
% the notion is more general.
% For instance, we could imagine reasoning about
% bias of a routing protocol, a spam filter, etc.
% In this section, we illustrate our verification technique and
% apply it to verifying group fairness of
% a simple hiring example.
% Despite the potential shortcomings of group fairness properties~\cite{dwork12},
% their simple formulation serves well as a proof of concept
% for our technique.

%\subsection{Illustrative Example}
\paragraph{A simple verification problem}
Consider the two programs in Figure~\ref{fig:ex}(a).
The program \texttt{popModel} is a probabilistic
program describing a simple model of the population.
% We view the population model as a \emph{probabilistic precondition},
% as it defines a joint probability distribution over the inputs of \texttt{dec},
% the decision-making program.
%
Here, a member of the population has three
attributes, all of which are real-valued:
\rone \texttt{ethnicity};
\rtwo \texttt{colRank}, the ranking of the college
the person attended (lower is better);
and \rthree \texttt{yExp}, the years of work experience
a person has.
%
We consider a person is a member of a protected group
if \texttt{ethnicity > 10}; we call this the
\emph{sensitive condition}.
%
The population model can be viewed as a \emph{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 \texttt{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 \emph{fairness target}).
The program implements a decision tree,
perhaps one generated by a machine-learning algorithm.
A person is hired if they attended a \emph{top-5}  college (\texttt{colRank <= 5}) or  have lots of experience compared to their college's ranking (\texttt{expRank > -5}).
Observe that \texttt{dec} \emph{does not access an applicant's ethnicity}.

Our goal is to establish
whether
the hiring algorithm \texttt{dec}  discriminates against
members of the protected minority.
%Concretely, to prove fairness,
Concretely, we attempt to prove the following property:
\[
   \frac{\pr{\texttt{hire}  \mid \texttt{min}}}
    {\pr{\texttt{hire} \mid \neg \texttt{min}}}
>
1 - \epsilon
\]
where \texttt{min} is shorthand for the sensitive condition
\texttt{ethnicity > 10},
and  $\epsilon$ is a small parameter set to $0.1$ for illustration.
% \footnote{Feldman et al.~\cite{feldman15} introduced this group fairness variation as an encoding of the \emph{Equality of Employment
% Opportunity Commission's}~\cite{EEOC} recommendation
% in the US.
% We could additionally impose a lower bound on the ratio
% to ensure that selection rate of the two groups is similar.
% }
%
Despite the potential shortcomings of this group fairness property~\cite{dwork12},
its simple formulation serves well as an illustration of our technique.

We can rewrite the above statement to eliminate
conditional probabilities as follows:
%
\begin{align}\label{al:fair}
\frac{\pr{\texttt{hire} \land \texttt{min}} \cdot \pr{\neg \texttt{min}}}
    {\pr{\texttt{hire} \land \neg \texttt{min}} \cdot \pr{\texttt{min}}}
>
1 - \epsilon
\end{align}
%
Therefore, to prove the above statement,
we need to compute a value for each of the probability
terms: $\pr{\texttt{hire} \land \texttt{min}}$, $\pr{\texttt{min}}$,
and $\pr{\texttt{hire} \land \neg \texttt{min}}$.
(Note that $\pr{\neg \texttt{min}} = 1 - \pr{\texttt{min}}$.)
%Similarly, to prove that the algorithm does \emph{not} satisfy the group fairness property,
%we attempt to prove
%that inequality~\ref{al:fair} does \emph{not} hold.

\begin{figure}[t]
  \centering
\fbox{\includegraphics[scale=1.0]{Ex}}
\caption{Simple illustrative example}
\label{fig:ex}
\end{figure}

Notice that, to prove or disprove inequality~\ref{al:fair},
all we need are good enough bounds on the values of such probabilities
and not their  exact values.
%Specifically, to prove the property, we need to compute a
%lower bound on the ratio
%$$\frac{\pr{\texttt{hire} \land \texttt{min}} \cdot \pr{\neg \texttt{min}}}
%    {\pr{\texttt{hire} \land \neg \texttt{min}} \cdot \pr{\texttt{min}}}$$
%If the lower bound is greater than $1-\epsilon$,
%then we know that the exact value of the ratio is also greater than $1-\epsilon$,
%and therefore the program satisfies group fairness.
%This can be performed by computing a lower bound for the numerator probabilities
%and an upper bound for the denominator probabilities.
%Conversely, to disprove the property, we need an upper bound on the ratio
%that is less than $1-\epsilon$.

For the purposes of illustration, we shall focus
our description on computing $\pr{\texttt{hire} \land \neg \texttt{min}}$.

\paragraph{Probabilistic verification conditions}
To compute the probability $\pr{\texttt{hire} \land \neg \texttt{min}}$,
we need to reason about the \emph{composition} of the two
programs, $\texttt{dec}\circ\texttt{popModel}$.
That is, we want to compute the probability that
\rone \texttt{popModel} generates a non-minority applicant,
and \rtwo \texttt{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 \emph{verification
condition} (\abr{VC}) generation for loop-free program fragments.

First, we encode \texttt{popModel} as follows:
\begin{align*}
\varphi_\texttt{pop} \equiv&~
\emph{ethnicity} > 10 \Rightarrow \emph{colRank}_1 = \emph{colRank} + 5 \\
&\land \emph{ethnicity} \leq 10 \Rightarrow \emph{colRank}_1 = \emph{colRank} \end{align*}
where subscripts are used to encode multiple occurrences
of the same variable (i.e., \abr{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 \texttt{dec} as follows (after simplification):
\begin{align*}
\varphi_\texttt{dec} \equiv&~
  \emph{expRank} = \emph{yExp}^i - \emph{colRank}^i\\
  &\land
  \emph{hire} \iff (\emph{colRank}^i \leq 5 \lor \emph{expRank} > -5)
\end{align*}
where variables with the superscript $i$ are the input
arguments to $\texttt{dec}$.
Now, to encode the composition $\texttt{dec} \circ \texttt{popModel}$,
we simply conjoin the two formulas---$\varphi_\texttt{pop}$ and $\varphi_\texttt{dec}$---and add equalities
between returns of \texttt{popModel} and arguments of \texttt{dec}.
$$\semprog \equiv \varphi_\texttt{pop} \land \varphi_\texttt{dec} \land
\emph{yExp}^i = \emph{yExp} \land \emph{colRank}^i = \emph{colRank}_1$$

Our goal is to compute the probability that a
non-minority applicant gets hired.
%Note that the probabilistically assigned variables
%are $\vars_p = \{\emph{ethnicity}, \emph{colRank},\emph{yExp}\}$.
%Therefore, we are asking,
Formally, we are asking,
\emph{what is the probability that the following formula is satisfied?}
$$\varphi \equiv \exists V_d \ldotp \semprog  \land \emph{hire} \land \emph{ethnicity}\leq10$$
 $V_d$ is the set of variables
that are not probabilistically assigned to, that is, all variables \emph{other than}
the three variables
$V_p = \{\emph{ethnicity}, \emph{colRank},\emph{yExp}\}$.
Intuitively, by \emph{projecting out} all non-probabilistic
variables, we get a formula $\varphi$ whose models are
the set of all probabilistic samplings
that lead to a non-minority applicant being generated and hired.

\paragraph{Weighted volume computation}
To compute the probability that $\varphi$
is satisfied, we begin by noting that
$\varphi$ is, geometrically,
a region in $\mathds{R}^3$, because it has three free real-valued variables, $V_p$.
The region $\varphi$ is partially illustrated in Figure~\ref{fig:ex}(b).
Informally, the probability of satisfying $\varphi$
is the probability of drawing values for the dimensions in $V_p$
that end up \emph{falling} in the region $\varphi$.
Therefore, the probability of satisfying $\varphi$
is its \emph{volume} in $\mathds{R}^3$,
\emph{weighted} by the probability density of each of the
three variables.
Formally:
%
$$\textstyle \pr{\texttt{hire} \land \neg \texttt{min}} = \int_{\varphi} p_ep_yp_c ~dV_p$$
%
where, e.g., $p_e$ is the \emph{probability density function}
of the distribution \texttt{gauss(0,10)}---the distribution
from which the value of \texttt{ethnicity} is drawn  in line 2 of \texttt{popModel}. Specifically, $p_e$ is a function of $\emph{ethnicity}$,
namely, $p_e(\emph{ethnicity}) = \frac{1}{10\sqrt{2\pi}}e^{-\frac{\emph{ethnicity}^2}{200}}$.


The primary difficulty here is that the region of integration is specified by an
arbitrary \abr{SMT} formula over an arithmetic theory.
So, \emph{how do we compute a numerical value for this integral?}
We make two interdependent observations:
\rone if the formula represents a \emph{hyperrectangular region}
in $\mathds{R}^n$---i.e., a box---then integration is typically simple,
due to the constant upper/lower bounds of all dimensions;
\rtwo we can \emph{symbolically decompose} an \abr{SMT}
formula into an (infinite) set of hyperrectangles.

Specifically, given our formula $\varphi$,
we  construct a new formula,
$\rect_\varphi$, where each model $m \models \rect_\varphi$
corresponds to a hyperrectangle that
underapproximates $\varphi$.
%is subsumed by $\varphi$, i.e., an \emph{underapproximation} of $\varphi$.
Therefore, by systematically finding disjoint hyperrectangles
inside of $\varphi$ and computing
their weighted volume, we iteratively improve a
\emph{lower bound} on the exact weighted volume of $\varphi$.
%One can view this process as a symbolic interpretation of
%\emph{Riemann integration} for arithmetic \abr{SMT} theories.
Figure~\ref{fig:ex}(c) shows a possible
underapproximation of $\varphi$ composed of four hyperrectangles.
Sec.~\ref{sec:vol} formalizes this technique  and proves
its convergence for decidable arithmetic theories.

%Figure~\ref{fig:ex}(c) shows a possible
%underapproximation of $\varphi$ composed of a number
%of hyperrectangles. Observe how  $\varphi$ is underapproximated by four hyperrectangles---%
%one large one and three hyperrectangles forming a ladder that underapproximates
%the slanted face of $\varphi$.
%
%Note that,
%if we are interested in an upper bound on the weighted volume of $\varphi$,
%we can simply compute a lower bound for the volume of the negation of  $\varphi$ and subtract it from 1. (Intuitively, the weighted volume of
%$\neg \varphi$ is  $1-\pr{\texttt{hire} \land \neg\texttt{min}}$.)

\paragraph{Proofs of group fairness}
We demonstrated how our technique reduces the problem
of computing probabilities to weighted volume computation.
%However, as we have shown, our technique computes
%a lower or an upper bound on the weighted volume of a given Boolean
%program property.
%For the purposes of proving fairness,
%we do not need an exact solution for the probabilities---all we need
%is a good enough solution that allows us to prove inequality~\ref{al:fair}.
%\loris{add disprove}
%Specifically, to prove fairness, we need to compute a
%lower bound on the ratio
%$$\frac{\pr{\texttt{hire} \land \texttt{min}} * \pr{\neg \texttt{min}}}
%    {\pr{\texttt{hire} \land \neg \texttt{min}} * \pr{\texttt{min}}}$$
%If the lower bound is greater than $1-\epsilon$,
%then we know that the exact value of the ratio is also greater than $1-\epsilon$,
%and therefore the program is fair.
%This can be performed by computing a lower bound for the numerator probabilities
%and an upper bound for the denominator probabilities.
%Conversely, to prove unfairness, we need an upper bound on the ratio
%that is less than $1-\epsilon$.
%
Figure~\ref{fig:ex}(d) illustrates a run of  our tool, \fs,  on this example.
\fs iteratively
improves lower and upper bounds for the
probabilities in the ratio,
and, therefore,
%it keeps refining a lower and an upper bound on the ratio.
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 $\epsilon = 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 \texttt{popModel}.
Looking at \texttt{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 \texttt{dec}
to \lstinline{expRank $\gets$ 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 \fs on the modified \texttt{dec} is illustrated
in Figure~\ref{fig:ex}(e), where the lower bound on the ratio
exceeds 0.9, thus proving our fairness property.

%!TEX root = paper.tex

impl.tex
\section{Implementation and Evaluation}\label{sec:impl}

%!TEX root = paper.tex

\newcommand{\dt}[1]{\abr{DT}$_{#1}$}
\newcommand{\dta}[1]{\abr{DT}$_{#1}^{\alpha}$}
\newcommand{\svm}[1]{\abr{SVM}$_{#1}$}
\newcommand{\svma}[1]{\abr{SVM}$_{#1}^{\alpha}$}
\newcommand{\nn}[2]{\abr{NN}$_{#1,#2}$}

% We implemented our algorithm in a tool called \fs.
% \fs was executed on a number of benchmarks, run
% on an Intel Core i5 3.20GHz CPU with 6 GB of RAM.

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 $\prog$ satisfies the condition
\[
    \frac{\pr{\prog(\vec{v}) = \true  \mid \emph{min}(\vec{v}) = \true \land \emph{qual}(\vec{v}) = \true}}
    {\pr{\prog(\vec{v}) = \true  \mid \emph{min}(\vec{v}) = \false \land \emph{qual}(\vec{v}) = \true}}
>
0.85
\]
We implement our algorithms in a tool called \fs
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 \fs to other probabilistic inference tools.

\subsection{Implementation}\label{ssec:impl}
We implemented our presented algorithms in a new tool
called \fs, which
employs  Z3~\cite{de2008z3} for \abr{SMT} solving and
Redlog~\cite{redlog} for quantifier elimination.
%In its current state, \fs supports step functions and Gaussians for probability distributions.
To compute the group fairness ratio, we decompose the conditional probabilities
into four joint probabilities:
$\pr{\prog(\vec{v}) \land \emph{min}(\vec{v}) \land \emph{qual}(\vec{v})}$,
$\pr{\prog(\vec{v}) \land \neg \emph{min}(\vec{v}) \land \emph{qual}(\vec{v})}$,
$\pr{\emph{min}(\vec{v}) \land \emph{qual}(\vec{v})}$, and
$\pr{\neg \emph{min}(\vec{v}) \land \emph{qual}(\vec{v})}$.
\fs 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.
%To compute the fairness
%ratio $\frac{\pr{\varphi_f \mid \varphi_s}}{\pr{\varphi_f \mid \lnot \varphi_s}}$,
%we use our symbolic volume computation method
%to compute lower bounds on the following six quantities:
%$\pr{\varphi_s}$,
%$\pr{\lnot \varphi_s}$,
%$\pr{\varphi_f \land \varphi_s}$,
%$\pr{\lnot(\varphi_f \land \varphi_s)}$,
%$\pr{\varphi_f \land \lnot \varphi_s}$, and
%$\pr{\lnot(\varphi_f \land \lnot \varphi_s)}$.
%These are sufficient to compute all quantities utilized
%in the algorithm shown in Figure~\ref{alg:fairsquare},
%including the upper bounds.
% for example,
% $\pru{\varphi_s} = 1 - \prl{\lnot \varphi_s}$.
% Note that we perform quantifier elimination
% on the formulas before computing the volume.
%
A \emph{round} of sampling involves \rone obtaining a
sample (hyperrectangle) for each of these eight quantities,
 \rtwo computing these samples' weighted volumes,
  \rthree
 updating the bounds on the group fairness ratio,
 and \rfour checking if bounds are precise enough
 to conclude \emph{fairness} or \emph{unfairness}.
Rounds of sampling are performed until a proof is found
or a timeout is reached.

% Recall this sample is Z3's solution of
% upper and lower bounds for all of the variables,
% forming a hyperrectangle.
% We use Scipy~\cite{scipy} to evaluate the integral
% of each variable with respect to its distribution
% over these bounds to obtain the weighted volume of the sample.
% After each round, the bounds on the fairness ratio
% are computed, and if they are precise enough
% to be conclusive (returning \emph{fair} or \emph{unfair}),
% the algorithm terminates.
% Otherwise, it arbitrarily continues with further rounds of sampling.

\paragraph{Sample maximization}
A key optimization is the maximization of hyperrectangles
obtained during sampling.
% All hyperrectangles obtained from the satisfiability queries
% have finite bounds.
% We check that hyperrectangle $H$ is contained in the
% formula $\varphi$ by ensuring the query
% $H \land \lnot \varphi$ is unsatisfiable;
% taking \abr{UNSAT} cores from this implication
% gives us see a set of potentially unconstrained bounds,
% which still preserve the implication $H \implies \varphi$.
% This set allows us to drop the unneccesary bounds, i.e.
% gives us hyperrectangles that extend infinitely.
%, as well as \abr{UNSAT} cores,
We use Z3's optimization capability
 to maximize and minimize the finite bounds
of all hyperrectangles, while still
satisfying the formula $\Psi$ (in Figures~\ref{alg:vol} and~\ref{fig:algadf}).
This process is performed greedily by extending a hyperrectangle in
one dimension at a time to find a
maximal hyperrectangle.
%\loris{it's not guaranteed to be maximal---Sam: it is \emph{maximal}, but it is not \emph{the maximum}}
If a dimension extends to infinity, then we drop that bound,
thus resulting in an unbounded hyperrectangle.
%The importance of sample maximization is illustrated in Figure~\ref{fig:opt}.
%\todo{too soon; we haven't discussed benchmarks yet}

\subsection{Benchmarks}
We trained a variety of machine-learning models on a popular income dataset~\cite{dataset}
used in related research on algorithmic fairness~\cite{feldman15,zemel13,calders10}
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:
\emph{high} ($>\$50{,}000$) or \emph{low}.
%The \emph{fairness target} is \emph{high salary} ($>\$50{,}000$)
%and the \emph{sensitive condition} is whether the \emph{sex} of
%the applicant is \emph{female}.
We would like to verify whether
salary decisions are fair to qualified female employees.
Using the Weka machine learning suite~\cite{weka},
we learned 11 different decision-making programs
(see, e.g., Bishop's textbook~\cite{bishop2006pattern}
for background)
which are listed in Figure~\ref{fig:summary}:

\begin{itemize}
\setlength{\itemsep}{-2pt}
\item Four \emph{decision trees}, named \dt{n},
    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.
\item Four \emph{support vector machines} with linear kernels, named \svm{n},
    where $n$ is the number of variables in the linear separator.
\item Three  \emph{neural networks} using \emph{rectified linear units}~\cite{nair2010rectified}, named \nn{n}{m},
    where $n$ is the number of input variables,
    and $m$ is the number of nodes in the single hidden layer.
\end{itemize}
%
% Note that \rone the support vector machines were quite sensitive to outliers
% in the data and were trained on a refined dataset that excluded outliers, and
% \rtwo the neural networks we prone to learning a model that
% did not hire for all inputs, so the coefficients were manually tweaked
% to make the programs more interesting.

As we will show in the next section,
some of these programs do not satisfy group fairness.
We introduced modifications of \dt{16} and \svm{4},
called \dta{16} and \svma{4},
that implement rudimentary forms of \emph{affirmative action}
for female applicants.
For \dta{16}, there is a 15\% chance it will flip a decision to give the low salary;
for \svma{4}, the linear separator is moved to increase the likelihood of hiring.

Additionally, we used three different probabilistic population models,
programs that define the probabilistic inputs, that were inferred from the same dataset:
\rone a set of \emph{independently distributed} variables,
\rtwo a \emph{Bayesian network} using a simple graph structure, and
\rthree the same Bayesian network,
but with an integrity constraint in the form of
an inequality between two of the variables.
%\begin{itemize}
%  \setlength{\itemsep}{-2pt}
%    \item A set of \emph{independently distributed} variables:
%        distributions of categorical variables are encoded as step functions,
%        while distributions of real-valued variables are approximated as Gaussians.
%    \item A \emph{clustered model}:
%        the dataset was split into two \emph{clusters} using the $k$-means algorithm.
%        Each cluster was fitted to a set of independent distributions, as
%        described above; the set of distributions to use is selected
%        by a coin toss weighted by the number of data points in each cluster.
%    \item A \emph{Bayesian network}:
%        we assumed a simple graph structure for the conditional dependencies
%        between variables (see Figure~\ref{fig:pop}).
%        The parameters of the network are independent
%        distributions learned from the dataset,
%        as in the above two cases.
%\end{itemize}
%
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
\iftoggle{long}{
(see Appendix~\ref{app:code} for an example).
}
{
(see \emph{supplementary materials} for an example).
}
%\begin{figure}[t]
%    \smaller
%    \begin{minipage}{0.1\textwidth}
%        \centering
%        \begin{tikzpicture}[every node/.style={draw}]
%            \path[use as bounding box, draw=none] (-0.375, -0.125) rectangle (0.875,-1.5);
%            \node(3) at (-0.125,-1.25) [circle] {};
%            \node(4) at (0.625,-1.25) [circle] {};
%            \draw[rounded corners=5pt,dashed] (-0.375, -1.5) rectangle (0.875, -1);
%            \node[draw=none] at (0.25,-1.25) {\ldots};
%        \end{tikzpicture}\\
%        Independent
%    \end{minipage}\hfill
%    \begin{minipage}{0.1\textwidth}
%        \centering
%        \begin{tikzpicture}[every node/.style={draw}]
%            \path[use as bounding box, draw=none] (-0.375, -0.125) rectangle (0.875,-1.5);
%            \node(1) at (0.25,-0.25) [circle] {};
%            \node[draw=none,anchor=east] at (0.125, -0.25) {(cluster)};
%            \node(3) at (-0.125,-1.25) [circle] {};
%            \node(4) at (0.625,-1.25) [circle] {};
%            \node(5)[draw=none] at (0.25,-1.125) {};
%            \draw[rounded corners=5pt,dashed] (-0.375, -1.5) rectangle (0.875, -1);
%            \node[draw=none] at (0.25,-1.25) {\ldots};
%            \draw [->,>=stealth] (1) edge (5);
%        \end{tikzpicture}\\
%        Clusters
%    \end{minipage}\hfill
%    \begin{minipage}{0.1\textwidth}
%        \centering
%        \begin{tikzpicture}[every node/.style={draw}]
%            \path[use as bounding box, draw=none] (-0.375, -0.125) rectangle (0.875,-1.5);
%            \node(1) at (0,-0.25) [circle] {};
%            \node[draw=none,anchor=east] at (-0.125, -0.25) {(sex)};
%            \node(2) at (0.5,-0.625) [circle] {};
%            \node[draw=none,anchor=west] at (0.625,-0.625) {(capital-gain)};
%            \node(3) at (-0.125,-1.25) [circle] {};
%            \node(4) at (0.625,-1.25) [circle] {};
%            \node(5)[draw=none] at (0,-1.125) {};
%            \node(6)[draw=none] at (0.5,-1.125) {};
%            \draw[rounded corners=5pt,dashed] (-0.375, -1.5) rectangle (0.875, -1);
%            \node[draw=none] at (0.25,-1.25) {\ldots};
%            \draw [->,>=stealth] (1) edge (2);
%            \draw [->,>=stealth] (1) edge (5);
%            \draw [->,>=stealth] (2) edge (6);
%        \end{tikzpicture}\\
%        Bayes Net
%    \end{minipage}\hspace{1cm}
%    \caption{Visualizations of \emph{conditional dependencies}
%    of variables (nodes) in our population models.}
%\label{fig:pop}
%\end{figure}


\subsection{Evaluation}

\begin{figure*}[t]
\centering
\begin{minipage}[b]{.6\textwidth}
\centering
\scriptsize
\setlength{\tabcolsep}{2.5pt}
\newcommand{\fair}[3]{\multicolumn{1}{c}{\cellcolor{green!20}\ding{51}} & { #1} & { #2} & { #3}}
\newcommand{\unfair}[3]{\multicolumn{1}{c}{\cellcolor{red!25}\ding{55}} & { #1} & { #2} & { #3}}
\newcommand{\stime}[4]{{\tiny $\overset{#1}{#2}$} & { #3} & { \abr{TO}} & { #4}}
\newcommand{\qtime}{\abr{TO}$_q$ & - & - & { \abr{TO}}}
\newcommand{\tblrow}[5]{#1 & #2 & #3 & #4 & #5}
\newcommand{\hhlinesep}{\hhline{|=##=##===##===##===|}}
\newcommand{\trih}{{ \it Res} & { \it \#} & { \it Vol} & { \it QE}}
\newcommand{\thline}{\specialrule{.15em}{.00em}{.00em}}
\renewcommand{\thline}{\hline}
{\def\arraystretch{1.0}
\begin{tabular}{cccccccccccccc}
    \toprule

    \multirow{3}{1.3cm}{\centering  \textbf{Decision program}} & \multirow{3}{*}{ \bf Acc} & \multicolumn{12}{c}{ \bf Population Model} \\
    & & \multicolumn{4}{c}{ Independent} & \multicolumn{4}{c}{ Bayes Net 1} &\multicolumn{4}{c}{ Bayes Net 2} \\
    \cmidrule(lr){3-6}
    \cmidrule(lr){7-10}
    \cmidrule(lr){11-14}
    & & \trih & \trih & \trih \\
    \midrule
    \tblrow{\dt{4}}{0.79}{\fair{10}{1.3}{0.5}}{\unfair{12}{2.2}{0.9}}{\unfair{18}{6.6}{2.2}} \\
    \tblrow{\dt{14}}{0.71}{\fair{20}{4.2}{1.4}}{\fair{38}{52.3}{11.4}}{\fair{73}{130.9}{33.6}} \\
    \tblrow{\dt{16}}{0.79}{\fair{21}{7.7}{2.0}}{\unfair{22}{15.3}{6.3}}{\unfair{22}{38.2}{14.3}} \\
    \tblrow{\dta{16}}{0.76}{\fair{18}{5.1}{3.0}}{\fair{34}{32.0}{8.2}}{\fair{40}{91.0}{19.4}} \\
    \tblrow{\dt{44}}{0.82}{\fair{55}{63.5}{9.8}}{\unfair{113}{178.9}{94.3}}{\unfair{406}{484.0}{222.4}} \\ \thline
    \tblrow{\svm{3}}{0.79}{\fair{10}{2.6}{0.6}}{\unfair{10}{3.7}{1.7}}{\unfair{10}{10.8}{6.2}} \\
    \tblrow{\svm{4}}{0.79}{\fair{10}{2.7}{0.8}}{\unfair{18}{13.3}{3.1}}{\unfair{14}{33.7}{20.1}} \\
    \tblrow{\svma{4}}{0.78}{\fair{10}{3.0}{0.8}}{\fair{22}{15.7}{3.2}}{\fair{14}{33.4}{63.2}} \\
    \tblrow{\svm{5}}{0.79}{\fair{10}{8.5}{1.3}}{\unfair{10}{12.2}{6.3}}{\qtime} \\
    \tblrow{\svm{6}}{0.79}{\stime{0.02}{35.3}{634}{2.4}}{\stime{0.09}{3.03}{434}{12.8}}{\qtime} \\ \thline
    \tblrow{\nn{2}{1}}{0.65}{\fair{78}{21.6}{0.8}}{\fair{466}{456.1}{3.4}}{\fair{154}{132.9}{7.2}}\\
    \tblrow{\nn{2}{2}}{0.67}{\fair{62}{27.8}{2.0}}{\fair{238}{236.5}{7.2}}{\fair{174}{233.5}{18.2}}\\
    \tblrow{\nn{3}{2}}{0.74}{\stime{0.03}{674.7}{442}{10.0}}{\stime{0.00}{5.24}{34}{55.9}}{\qtime} \\ \thline
\end{tabular}}
\end{minipage}%
\begin{minipage}[b]{.39\textwidth}
\centering
\pgfplotstableread{
    X   M    Tool DT  SVM NN
    1   Ind  F2   5   4   2
    2   BN1  F2   5   4   2
    3   BN2  F2   5   3   2
    5   Ind  \abr{PSI}  5   5   0
    6   BN1  \abr{PSI}  3   0   0
    7   BN2  \abr{PSI}  0   0   0
    9   Ind  \abr{VC}   4   0   0
    10  BN1  \abr{VC}   3   0   0
    11  BN2  \abr{VC}   0   0   0
}\datatable
\begin{tikzpicture}[baseline=(current bounding box.center)]
  \footnotesize
    \smaller
    \begin{axis}[
            width=6.5cm,
            xtick=data,
            xticklabels from table={\datatable}{M},
            xticklabel style={rotate=90,anchor=east},
            xtick style={draw=none},
            enlarge x limits=0.1,
            ymin=0,ymax=13,
            ytick={0,2,4,6,8,10,12},
            ybar stacked,
            bar width=7pt,
            legend style={draw=none,legend columns=-1},
            after end axis/.append code={
                \path
                (axis cs:2,0) node [yshift=-7ex] {\fs}
                (axis cs:6,0) node [yshift=-7ex] {\abr{PSI}~\cite{Gehr16}}
                (axis cs:10,0) node [yshift=-7ex] {\abr{VC}~\cite{sankaranarayanan13}};
            }
        ]
        \addplot[draw=black,fill=black!20] table[x=X,y=DT]\datatable; \addlegendentry{DT}
        \addplot[draw=black,fill=black!40] table[x=X,y=SVM]\datatable; \addlegendentry{SVM}
        \addplot[draw=black,fill=black!60] table[x=X,y=NN]\datatable; \addlegendentry{NN}
    \end{axis}
\end{tikzpicture}
\end{minipage}
\caption{
{\bf(Left)}
Results of \fs applied to 39 group fairness problems.
\emph{Res} indicates
\ding{51} for fair; \ding{55} for unfair.
\emph{Vol} indicates the time (s)
of the sampling procedure;
\# is the number of \abr{SMT} calls.
If sampling timed out (900s),
\emph{Res} denotes the latest bounds
on the fairness ratio.
\emph{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 \abr{TO}$_q$ for \emph{Res}.
\emph{Acc} is the training set accuracy
for each of the programs.
{\bf(Right)}
Comparison of the number of benchmarks
that \fs, \abr{PSI}~\cite{Gehr16}, and \abr{\abr{VC}}~\cite{sankaranarayanan13} were able to
solve.
}
\label{fig:summary}
\end{figure*}

%\begin{figure}[t]
%    \centering
%    \small
%    \begin{minipage}{0.15\textwidth}
%        \centering
%        \includegraphics[scale=0.2]{quanfigs/exact.pdf} \\
%        (a)
%    \end{minipage}\hfill
%    \begin{minipage}{0.15\textwidth}
%        \centering
%        \includegraphics[scale=0.2]{quanfigs/asymp.pdf} \\
%        (b)
%    \end{minipage}\hfill
%     \begin{minipage}{0.15\textwidth}
%        \centering
%        \includegraphics[scale=0.2]{quanfigs/non.pdf} \\
%        (c)
%    \end{minipage}
%\caption{
%    Upper (red) and lower (blue) bounds of
%    fairness ratio vs \# of sampling rounds:
%    (a)  \dt{16},
%    (b)  \svm{3},
%    and (c) \svm{5}.
%    All under the clustered model.
%}
%\label{fig:behavior}
%\end{figure}

In this section, we discuss the ability of \fs
to verify qualified group fairness or unfairness for 39 problems,
as summarized in Figure~\ref{fig:summary}.
In these problems, $\emph{min}(\vec{v}) = \true$ when the applicant is female,
and $\emph{qual}(\vec{v}) = \true$ in two different scenarios:
first, we consider the case when $\emph{qual}$ is tautologically true,
and second, when the applicant is at least 18 years of age.
Figure~\ref{fig:summary} shows only the former case,
as while the numbers are different for the two cases,
the qualitative results are quite similar.
\iftoggle{long}{
See Appendix~\ref{app:res} for the full table of results.
}
{
See \emph{supplementary materials} for the full table of results.
}
We fix $\epsilon = 0.15$.
To guide volume computation,
all Gaussian distributions with mean $\mu$ and variance $\sigma^2$
use \abrs{ADF} with 5 equal-width steps spanning $(\mu - 3\sigma^2,\mu + 3\sigma^2)$---analogous to Figure~\ref{fig:approx}(a).
\fs was able to solve 32 of the 39 problems, proving 21 fair and 11 unfair,
as shown in the left table of Figure~\ref{fig:summary}.

%\todo{probably this paragraph not needed}
%Decision trees divide the input
%space into rectilinear boundaries, which are exactly a union
%of hyperrectangles.  Consequently, they are the best case for
%our tool, since we can perform exact volume computation.
%%
%Support vector machines, however, separate the space with
%an angled hyperplane.
%As illustrated in Figure~\ref{fig:decomp},
%potentially many samples are needed to obtain
%good approximations of the weighted volumes.
%(The number of hyperrectangles needed to approximate
%a region with an angular boundary
%also scales in the number of dimensions,
%i.e., in the number of variables.)
%Especially when the separating hyperplane intersects
%with regions of high probability,
%our tool requires a lot of samples and thus
%takes longer to compute fairness for these models.
%%
%The rectified linear neural networks are essentially
%a collection of support vector machines integrated
%into a decision tree, so it is expected that
%computing fairness for them is generally harder
%than the support vector machine instances.

Consider the results for \dt{4}:
\fs 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 \abr{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 \dt{44}
under the Bayes Net 1 population model:
\fs was unable to conclude fairness or unfairness
after 900 seconds of volume computation
(denoted by \abr{TO} in the \emph{Vol} column).
The lower and upper bounds of the fairness ratio
it had computed at that time are listed in the \emph{Res} column:
in this case, the value of the fairness ratio
is within $[0.70,0.88]$, which is not
precise enough for the $\epsilon = 0.15$ requirement
(but would be precise enough for $\epsilon$
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 \emph{sex},
yet several are unfair.
%This result is supported by existing work
%that shows a \emph{na\"{i}ve Bayes} classifier for this
%dataset has substantial bias
%with respect to sex~\cite{feldman15}.

Figure~\ref{fig:summary} also shows that the modifications
for affirmative action in \dta{16} and \svma{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, \fs is powerful enough to reason about
group fairness for many non-trivial programs.

%Figure~\ref{fig:behavior} illustrates three patterns
%of convergence for fairness ratio
%computation in \fs.  The cases are:
%  (a) the exact value of the ratio is eventually computed;
%    (b) the upper and lower bounds of the computation quickly
%        become asymptotically close to the exact value; and
%    (c) the upper and lower bounds of the computation
%        converge very slowly to the exact value.
%%
%Recall that it is sufficient
%to prove merely that the upper or lower bound of the fairness ratio
%falls below or above, respectively, $1 - \epsilon$.
%Consequently, the result for case (c) can be conclusive,
%but is sensitive to the given $\epsilon$.

\subsection{Comparison to other tools}
We ran our benchmarks on the two other \emph{exact} probabilistic inference tools
that accept the same class of problems,
and report these results in the right bar graph in Figure~\ref{fig:summary}.
First, we compare to Sankaranarayanan et al.'s tool~\cite{sankaranarayanan13} (\abr{VC}),
which 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 \abr{PSI}~\cite{Gehr16},
which symbolically computes representations of
the posterior distributions of variables.
When closed form solutions to the \abrs{CDF} exist,
\abr{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 \abr{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 \abr{VC}
were often not tight enough for a proof.
%It should be noted that when these tools did terminate,
%it was often on the order of seconds;
\iftoggle{long}{
See Appendix~\ref{app:res} for the full quantitative table of results.
}
{
See \emph{supplementary materials} 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 \emph{constant}.
However, the presence of inequalities
involving multiple variables
can result in two phenomena:
\rone the lack of closed form posterior \abrs{CDF},
as reflected in the output of \abr{PSI}, and
\rtwo angled boundaries in the decision space
that are hard to approximate with hyperrectangles.
These inequalities occur in the \abrs{SVM},
neural networks, and the Bayes Net 2 population model.
Consequently,
\abr{VC} fails to produce good bounds in these cases;
\fs performs better, but also fails on the largest \abr{SVM}.

\subsection{Effect of Parameters}

\begin{figure}[t]
  \smaller
  \centering
    \setlength{\tabcolsep}{0pt}
    \begin{tabular}{rcc}
        & \dt{16} & \svm{4} \\
        \rotatebox{90}{\smaller \hspace{1em} \parbox{3cm}{\centering fairness ratio without sample maximization}} &
        \begin{overpic}[scale=0.25]{quanfigs/DT_nomax_dist.pdf} \put(80,50) {(a)} \end{overpic}&
        \begin{overpic}[scale=0.25]{quanfigs/SVM_nomax_dist.pdf} \put(80,50) {(c)} \end{overpic} \\
        \rotatebox{90}{\smaller \hspace{1.5em} \parbox{3cm}{\centering fairness ratio with sample maximization}} &
        \begin{overpic}[scale=0.25]{quanfigs/DT_max_dist.pdf} \put(80,50) {(b)} \end{overpic}&
        \begin{overpic}[scale=0.25]{quanfigs/SVM_max_dist.pdf} \put(80,50) {(d)} \end{overpic}
    \end{tabular}
\caption{
    Fairness ratio vs. rounds of sampling
    for \dt{16} and \svm{4}
    (independent pop model)
    differing on \abrs{adf}
    and sample maximization.
    In (b) two runs end at the exact value.
    Outside the visible range are:
    (a)(c) upper and lower bounds of \emph{uniform} and \emph{none};
    (d) upper bounds of \emph{none}.
}
\label{fig:opt}
\end{figure}

\begin{figure}[t]
    \centering
    \pgfplotsset{every tick label/.append style={font=\tiny}}
    \setlength{\tabcolsep}{0pt}
    \smaller
    \begin{tabular}{rcc}
        & ~~~~~~~~\dt{16} & ~~~~~~~~\svm{4} \vspace{0em} \\
        \rotatebox{90}{\smaller \hspace{1em} \parbox{2cm}{\centering average weighted volume per sample}} &
        \multicolumn{1}{r}{\begin{tikzpicture}
            \begin{axis}[
                    width=5cm,
                    height=3.2cm,
                    legend style={font=\tiny, legend columns=-1, at={(1,1.5)}},
                    ybar, ymin=0, ymax=0.15,
                    yticklabel style={
                        /pgf/number format/fixed,
                        /pgf/number format/precision=2
                    },
                    scaled y ticks=false,
                    symbolic x coords={5step,uniform,none},
                    xtick=data,
                    xticklabels={
                        5-step,
                        uniform,
                        none,
                    },
                    xticklabel style={text width=1cm, align=center},
                    bar width=7, x=1cm, enlarge x limits={abs=0.5cm},
                    after end axis/.append code={
                        \path (rel axis cs:0.9,0.85) node {(a)};
                    }
                    ]
                \addplot[ybar,fill=black] coordinates {
                    (5step, 0.1429)
                    (uniform, 0.1429)
                    (none, 0.0307)
                };
                \addplot[ybar,fill=black!20] coordinates {
                    (5step, 0.0236)
                    (uniform, 0.0123)
                    (none, 0.0022)
                };
                \legend{\dt{16} max, \dt{16} nomax}
            \end{axis}
        \end{tikzpicture}} &
            \multicolumn{1}{r}{\begin{tikzpicture}
            \begin{axis}[
                    width=5cm,
                    height=3.2cm,
                    legend style={font=\tiny, legend columns=-1, at={(1,1.5)}},
                    ybar, ymin=0, ymax=0.15,
                    symbolic x coords={5step,uniform,none},
                    yticklabel style={
                        /pgf/number format/fixed,
                        /pgf/number format/precision=2
                    },
                    scaled y ticks=false,
                    xtick=data,
                    xticklabels={
                        5-step,
                        uniform,
                        none,
                    },
                    xticklabel style={text width=1cm, align=center},
                    bar width=7, x=1cm, enlarge x limits={abs=0.5cm},
                    after end axis/.append code={
                        \path (rel axis cs:0.9,0.85) node {(c)};
                    }
                    ]
                \addplot[ybar,fill=black] coordinates {
                    (5step, 0.0331)
                    (uniform, 0.0330)
                    (none, 0.0206)
                };
                \addplot[ybar,fill=black!20] coordinates {
                    (5step, 0.0326)
                    (uniform, 0.0218)
                    (none, 0.0111)
                };
                \legend{\svm{4} max, \svm{4} nomax}
            \end{axis}
        \end{tikzpicture}}  \\
        \rotatebox{90}{\smaller \hspace{1.5em} \parbox{2cm}{\centering average time per round}} &
        \multicolumn{1}{r}{\begin{tikzpicture}
            \begin{axis}[
                    width=5cm,
                    height=3.2cm,
                    ybar, ymin=0, ymax=2,
                    symbolic x coords={5step,uniform,none},
                    xtick=data,
                    xticklabels={
                        5-step,
                        uniform,
                        none,
                    },
                    xticklabel style={text width=1cm, align=center},
                    bar width=7, x=1cm, enlarge x limits={abs=0.5cm},
                    after end axis/.append code={
                        \path (rel axis cs:0.9,0.85) node {(b)};
                    }
                    ]
                \addplot[ybar,fill=black] coordinates {
                    (5step, 0.9721)
                    (uniform, 0.8798)
                    (none, 0.6270)
                };
                \addplot[ybar,fill=black!20] coordinates {
                    (5step, 0.2491)
                    (uniform, 0.1924)
                    (none, 0.1908)
                };
            \end{axis}
        \end{tikzpicture}} &
        \multicolumn{1}{r}{\begin{tikzpicture}
            \begin{axis}[
                    width=5cm,
                    height=3.2cm,
                    ybar, ymin=0, ymax=2,
                    symbolic x coords={5step,uniform,none},
                    xtick=data,
                    xticklabels={
                        5-step,
                        uniform,
                        none,
                    },
                    xticklabel style={text width=1cm, align=center},
                    bar width=7, x=1cm, enlarge x limits={abs=0.5cm},
                    after end axis/.append code={
                        \path (rel axis cs:0.9,0.85) node {(d)};
                    }
                    ]
                \addplot[ybar,fill=black] coordinates {
                    (5step, 1.4835)
                    (uniform, 1.6445)
                    (none, 1.3091)
                };
                \addplot[ybar,fill=black!20] coordinates {
                    (5step, 0.3002)
                    (uniform, 0.2006)
                    (none, 0.2289)
                };
            \end{axis}
        \end{tikzpicture}} \\
    \end{tabular}
\caption{
    Effect of optimizations
    on  performance of \fs for
    \dt{16} and \svm{4} (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.
}
\label{fig:optbar}
\end{figure}

We now explore  the effects
of the approximate density functions (see Section~\ref{sec:density})
and of the sample maximization optimizations.
These results are captured in Figures~\ref{fig:opt} and~\ref{fig:optbar}.

There are three instances of \abrs{ADF} in Figure~\ref{fig:opt}
used to guide the sampling to high-probability regions:
\rone \emph{none} indicates that no \abr{ADF} is used, i.e.,
we used \volalg instead of \volopt;
\rtwo \emph{uniform} indicates that each $\texttt{gauss}(\mu, \sigma^2)$ is approximated by a
uniform function spanning $(\mu - 3\sigma^2, \mu + 3\sigma^2)$
(similar to Figure~\ref{fig:approx}(c)); and
\rthree \emph{5-step} indicates that each Gaussian is approximated by
a step function of 5 equal-width regions spanning that same domain
(similar to Figure~\ref{fig:approx}(a)).
%
Another variable, \emph{max} or \emph{nomax},
denotes whether the sample maximization optimization is enabled
(as described in Section~\ref{ssec:impl}).

Each combination of these techniques is run on two of our benchmarks:
\dt{16} and \svm{4} under the independent population model.
Figure~\ref{fig:opt}(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 \emph{uniform} and \emph{none}
are not even visible, as the bounds never fall within $[0.01,4.0]$.
Plots (b) and (d) show that when sample maximization \emph{is} employed,
the choice between the uniform and 5-step approximations is
not as substantial on these benchmarks, although
\rone the better approximation gets better bounds faster, and
\rtwo using none results in substantially worse bounds.

Figure~\ref{fig:optbar} plot (a) and (c) show that employing \abrs{ADF}
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 \abrs{ADF} and
sample maximization are both
 necessary for adequate performance of \fs.

intro.tex
\section{Introduction}\label{sec:intro}

% Probabilistic programs allow
% us to encode uncertainty through probabilistic assignments.
A number of very interesting applications of program analysis
have been explored in the probabilistic setting:
reasoning about cyber-physical systems~\cite{sankaranarayanan13},
proving differential privacy of complex algorithms~\cite{barthe14},
reasoning about approximate programs and hardware~\cite{sampson2014expressing,carbin13},
synthesizing control programs~\cite{chaudhuri14},
amongst many others.
In this paper, we turn our attention to the problem
of \emph{quantifying program bias}.

\paragraph{Program bias}
% Intuitively, \emph{program bias} tries to quantify \emph{how much}
% the output of the program  depends on particular
% subsets of values of its variables.
%Of course, we expect every program's output
%to be dependent on its variables.
Programs that make decisions can be biased.
Consider, for instance,
automatic grading of
writing prompts for standardized tests~\cite{attali06};
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~\cite{articleTimesHiring,articleWired},
welfare allocation~\cite{articleSlateWelfare},
prison sentencing~\cite{articlePropublica},
policing~\cite{articleGuardianCrime,perry2013predictive},
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 \emph{algorithmic fairness}
% which is ultimately
% a consequence of program bias,
has recently captured
the attention of a broad spectrum of
experts, within computer science and without~\cite{dwork12,zemel13,feldman15,calders10,datta2015automated,articlePropublica,articWsjstaples,sweeney2013discrimination,tutt2016fda,ajunwa2016hiring,barocas2014big}.

Fairness and justice
have always been a ripe topic for philosophical debate~\cite{rawls2009theory},
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~\cite{feldman15,ruggieri2014using,dwork12,hardt16,friedler16}.
Notable formulations of fairness include
\emph{individual fairness}, which dictates
that \emph{similar} inputs must result in \emph{similar} outputs,
and \emph{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 \emph{quantitative properties} of decision-making
 programs.

%---however, the fairness community
%has not yet come to a consensus on what the \emph{best} notion is.
%\todo{Citations?}

\paragraph{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.
% features of individuals,
% which are effectively random variables drawn from
% a joint distribution that represents the demographics of some population.
Quantifying program bias becomes a matter of reasoning about
probabilities of program executions.

Consider a hiring program $\prog$ that takes as input a vector
of arguments $\vec{v}$ representing a job applicant's record.
One of the arguments $v_s$ in the vector $\vec{v}$ states whether the person
is a member of a protected minority or not,
and similarly $v_q$ in $\vec{v}$ states whether
the person is qualified or not.
Evaluating $\prog(\vec{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 \emph{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:
%
\[
\frac{\pr{\prog(\vec{v}) = \true \mid v_s = \true \land v_q = \true}}
    {\pr{\prog(\vec{v}) = \true \mid v_s = \false \land v_q = \true}}
>
1 - \epsilon
\]
Here, $\epsilon$ is a small prespecified value.
In other words, the probability of hiring
a person $\vec{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 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 poses a number of challenges.
% Sampling-based methods to estimate  probabilities %in the formula
% produce bounds on the value that hold with high probability, which
% may not be suitable in, for example, a legal setting;
% furthermore, these methods cannot infer
% absolute bounds on conditional probabilities and  ratios of probabilities.
% Deductive techniques that can compute exact bounds on probabilities
% suffer from incompleteness~\cite{PSI}. %, which poses practical limitations.
%but were  not sufficient
%to prove our program bias assertions
%(see Section~\ref{sec:impl});

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 \emph{weighted volume} of
the logical encoding of a program in real arithmetic.
We then utilize a novel \emph{symbolic volume computation} algorithm
that exploits the power of \abr{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 \abr{SMT} with this expressivity and guarantees.
We implement our algorithm in a tool called \fs, which we
%an automated verification tool
%for quantifying program bias
%in terms of probabilities of outcomes and
evaluate on a number
  decision-making
  programs generated by a range of machine-learning
  algorithms from real-world data.

\paragraph{Contributions}
This paper makes a number of conceptual, algorithmic, and practical contributions:
\begin{itemize}
  \item
  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 \emph{real closed fields},
  that utilizes an \abr{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 \abr{SMT} with this generality and guarantees.
  (Sec.~\ref{sec:vol})

  \item We present an automated verification tool, \fs, 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.~\ref{sec:impl})

\end{itemize}

\iftoggle{long}{}{
The \emph{supplementary materials} include a full version
of this paper with proofs, as well as detailed experimental results
and the benchmark programs from Section~\ref{sec:impl}.
}

%!TEX root = paper.tex

paper.tex
\documentclass[10pt,nocopyrightspace,onecolumn,preprint,numbers]{sigplanconf}
\input{preamble.tex}
\toggletrue{long}

\begin{document}

%\setlength{\abovedisplayskip}{3pt}
%\setlength{\belowdisplayskip}{3pt}

\title{Quantifying Program Bias
}

\authorinfo{
Aws Albarghouthi, Loris D'Antoni, and Samuel Drews
}
           {
University of Wisconsin--Madison}
           {
           }


\authorinfo{ Aditya Nori
}
           {

Microsoft Research}
           {
           }


\maketitle
\begin{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 \emph{biased}.
  %When we automate decision-making tasks,
  %we risk instilling \emph{bias} in programs.
  We propose a
  novel probabilistic program analysis technique
  and apply it to quantifying bias in decision-making programs.
  Specifically, we \rone present a sound and complete
  automated verification technique for proving quantitative
  properties of probabilistic programs;
  \rtwo show that certain notions of bias,
  recently proposed in the \emph{fairness} literature,
  can be phrased as quantitative correctness properties; and
  \rthree present \fs, the first
  verification tool
  for quantifying program bias, and evaluate it on
  a range of
  decision-making
  programs.
  %  generated by a range of machine-learning
  % algorithms from real-world data.
% Algorithms have become powerful
% arbitrators of an ever-expanding range of sensitive
% decisions with significant societal
% impact---making hiring decisions,
% deciding credit scores, deciding welfare eligibility,
% predicting crime activity, directing policing,
% influencing public policy,
% amongst many others.
% And the indicators are all clear: Algorithmic decision-making
% is unlikely to stop in the foreseeable future;
% in fact, we are merely witnessing its very beginning.
%
% With the range and sensitivity of algorithmic decisions
% expanding at a break-neck speed,
% it is imperative that we aggressively investigate
%  \emph{algorithmic fairness
% and bias}.
% In this paper, we ask the following question:
% \emph{Does a decision-making program discriminate
% against a protected group?}
% Specifically,
% we \rone formalize algorithmic fairness as a
% probabilistic--program--verification problem,
% and \rtwo present a novel sound and complete
% automated verification technique, based on
% volume computation, for
% proving  (un)fairness of decision-making programs.
%
% We implemented our techniques in a new fairness
% verification tool called \fs, and applied
% it to proving fairness of decision-making
% programs generated by a range of machine learning
% algorithms from real-world data.
% Our results highlight the power and promise
% of our approach at attacking this new and
% increasingly important problem.
\end{abstract}


\input{intro}
\input{example}
\input{prelims}
\input{vol}
%\input{alg}
\input{impl}
%\input{discussion}
\input{related}
%\input{conclusion}

\bibliographystyle{abbrvnat}
\bibliography{biblio}

\appendix
\iftoggle{long}{
\input{proofs}
\input{code}
\input{results}
}

\end{document}

preamble.tex
\usepackage{etex}
\usepackage[usenames,dvipsnames,svgnames,table]{xcolor}
\usepackage{prooftree}
\usepackage{listings}
\usepackage{pifont}
\usepackage{multirow}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{xspace}
\usepackage{graphicx}
\usepackage{lastpage}
\usepackage{paralist}
\usepackage{fancyhdr}
\usepackage{wrapfig}
\usepackage{balance}
\usepackage{tikz}
\usepackage{ctable} % for \specialrule command
\usepackage{pgfplots}
\usepackage{amsthm}
\usepackage{dsfont}
\usepackage{colonequals}
\usepackage{enumitem}
\usepackage{wasysym}
\usepackage{mathtools}
\usepackage{listings}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{pifont}
\usepackage{stmaryrd}
\usepackage{comment}
\usepackage[hyphens]{url}
\usepackage{etoolbox}
\usepackage{pbox}
\usepackage{relsize}
\usepackage{csquotes}
\usepackage{scalerel}
\usepackage{setspace}
\usepackage[T1]{fontenc}
%\usepackage{euler}
\usepackage{manfnt}
\usepackage{tabularx}
\usepackage{hhline}
\usepackage{titlesec}
\usepackage{enumitem}% http://ctan.org/pkg/enumitem
\usepackage[percent]{overpic}

\setlength{\skip\footins}{6pt}


\titlespacing*{\section}
{0pt}{1.3ex plus .1ex minus .1ex}{1.ex plus .1ex}
\titlespacing*{\subsection}
{0pt}{1.3ex plus .1ex minus .1ex}{1.ex plus .1ex}

%\usepackage{mathabx}

%% code from mathabx.sty and mathabx.dcl
% getting widehat from mathabx, because importing
% the package changes math symbols
\DeclareFontFamily{U}{mathx}{\hyphenchar\font45}
\DeclareFontShape{U}{mathx}{m}{n}{
      <5> <6> <7> <8> <9> <10>
      <10.95> <12> <14.4> <17.28> <20.74> <24.88>
      mathx10
      }{}
\DeclareSymbolFont{mathx}{U}{mathx}{m}{n}
\DeclareFontSubstitution{U}{mathx}{m}{n}
\DeclareMathAccent{\widecheck}{0}{mathx}{"71}
\DeclareMathAccent{\wideparen}{0}{mathx}{"75}

\def\cs#1{\texttt{\char`\\#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%

\DeclareMathOperator*{\lp}{\scalerel*{\llparenthesis}{\sum}}
\DeclareMathOperator*{\rp}{\scalerel*{\rrparenthesis}{\sum}}


\newcommand{\ctx}{{\mbox{\scriptsize$\Gamma$}}}
\newcommand{\mkey}{\tau^\ctx}
\newcommand{\iffe}{\emph{iff}\xspace}
%short vs long paper toggle
\newtoggle{long}

%itemize
\def\labelitemi{--}



%\renewcommand{\Gamma}{\scalebox{.5}{\Gamma}}

%enumeration - description options

%theorems if needed
\theoremstyle{definition}
\newtheorem{example}{Example}
\AtEndEnvironment{example}{\null\hfill\qedsymbol}%
\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}
\newtheorem{task}{Task}
%\theoremstyle{definition}
\newtheorem{definition}{Definition}
%\AtEndEnvironment{definition}{\null\hfill\qedsymbol}%

\algrenewcommand\alglinenumber[1]{\smaller #1:}
% algorithm listings
\makeatletter
\lst@InstallKeywords k{attributes}{attributestyle}\slshape{attributestyle}{}ld
\makeatother
\newcommand\gtt[1]{{\color{gray}\texttt{\textbf{#1}}}}
\newcommand\nf[1]{{\normalfont#1}}
\usepackage[scaled=0.80]{DejaVuSansMono}
\lstdefinestyle{customc}{
  belowcaptionskip=1\baselineskip,
  breaklines=true,
  xleftmargin=\parindent,
  language=C,
  showstringspaces=false,
  basicstyle=\small\ttfamily,
  keywordstyle=\small\bfseries\ttfamily\color{NavyBlue},
  commentstyle=\itshape\color{black},
  identifierstyle=\ttfamily\color{black},
  stringstyle=\itshape\color{NavyBlue},
  keywords={ map, flatMap, reduce, return, define, gaussian, gauss, step, elif, then, skip, if, else, reduceByKey, sensitiveAttribute, fairnessTarget, qualified, def},
moreattributes={ let, where}, % etc...
attributestyle = \itshape\ttfamily\color{attributecolor}
}

\setlist[itemize]{
    %labelindent=.4cm,
    %style=unboxed,
    %leftmargin=.4cm,
    %font=\itshape,
    topsep=.5ex,
    itemsep=0ex,
    leftmargin=1em,
    %itemindent=1em
}

\setlist[description]{
    labelindent=.4cm,
    style=unboxed,
    leftmargin=.4cm,
    font=\itshape,
    topsep=.5ex,
    itemsep=0ex
}

\lstset{columns=fullflexible,
        mathescape=true,
        literate=,
        numbersep=5pt,
        numberstyle=\tt\color{gray}
}
\lstset{escapechar=@,style=customc}

%space above and below
\lstset{aboveskip=2pt,belowskip=2pt}

%\setlist[itemize]{topsep=3pt,itemsep=-2pt}
%\setlist[enumerate]{topsep=3pt,itemsep=0pt}
%\setlist[description]{
%    labelindent=.3cm,
%   style=unboxed,
%    leftmargin=.3cm,
%    font=\itshape,
%    topsep=3pt,
%    itemsep=-.5ex
%}


%change \paragraph looks and spacing
\renewcommand{\paragraph}[1]{\vspace{.03in}\noindent\textbf{\emph{#1}}~}

%subsubsection that is emph
\newcommand{\psec}[1]{\\\noindent\*{\text{\emph{#1}}}\\\noindent}


%footnote change to *
%\renewcommand*{\thefootnote}{\fnsymbol{footnote}}
%\setlength{\bibspacing}{\baselineskip}


%comments
\newcommand{\aws}[1]{{\color{RedViolet}{\noindent\textsf{\textbf{Aws}---#1}}}}
\newcommand{\an}[1]{{\color{Blue}{\noindent\textsf{\textbf{AN}---#1}}}}
\newcommand{\calvin}[1]{{\color{MidnightBlue}{\noindent\textsf{\textbf{Calvin}---#1}}}}
\newcommand{\todo}[1]{{\color{BrickRed}{\noindent\textsf{\textbf{{\scriptsize~TODO}}---#1}}}}
\newcommand{\loris}[1]{{\color{Blue}{\noindent\textsf{\textbf{{\scriptsize~LORIS}}---#1}}}}
\newcommand{\note}[1]{
{\color{darkgray}{\noindent{\it\sffamily \textbf{Notes:} #1}}}
}
%\renewcommand{\loris}[1]{}
%\renewcommand{\todo}[1]{}
%\renewcommand{\aws}[1]{}
%\renewcommand{\an}[1]{}


%page style
\newcommand{\shortenpar}{% try to make current paragraph one line shorter
  \looseness=-1%
}
\everypar{\looseness=-1}

%gray color
\definecolor{lightgray}{gray}{0.9}
\definecolor{midgray}{gray}{0.65}
\definecolor{darkgray}{gray}{0.3}

%dark blue

%gray box
%gray box
\newcommand\smallbox[1]{%
  {
  \begin{center}
  \centering
  \vskip\baselineskip%
  \par\noindent{%
    \begin{minipage}{.2\textwidth}#1\end{minipage}%
  }%
  \vskip\baselineskip%
  \end{center}
  }
}



\newcommand\graybox[1]{%
  {
  \vskip\baselineskip%
  \par\noindent\colorbox{lightgray}{%
    \begin{minipage}{.45\textwidth}
        \textbf{Notes:}
        #1
    \end{minipage}%
  }%
  \vskip\baselineskip%
  }
}

%use for abbreviations, e.g., UFO --> \abr{UFO}
\newcommand{\abr}[1]{\textsc{\MakeLowercase{#1}}}

\newcommand{\abrs}[1]{\abr{#1}{\footnotesize{s}}\xspace}

% vectors bold - arrows add a lot of clutter
\renewcommand{\vec}[1]{\boldsymbol{#1}}

% leq / geq with lower bar parallel to >
\renewcommand{\leq}{\leqslant}
\renewcommand{\geq}{\geqslant}


% dot for quantifiers
\newcommand*{\anddot}{%
  \mathclose{}%
  \nonscript\mskip1\thinmuskip
  \boldsymbol{.}%
  \;%
  \mathopen{}%
}

%qed is solid black as opposed to empty square
\renewcommand{\qedsymbol}{{\footnotesize{$\blacksquare$}}}

%in-text enumeration
\newcommand{\rone}{(\emph{i})~}
\newcommand{\rtwo}{(\emph{ii})~}
\newcommand{\rthree}{(\emph{iii})~}
\newcommand{\rfour}{(\emph{iv})~}
\newcommand{\rfive}{(\emph{v})~}

\newcommand{\smallparagraph}[1]{\noindent \emph{#1}}

%keyword color
\definecolor{keywordcolor}{gray}{0.0}
\definecolor{attributecolor}{gray}{0.0}
%wildcard color
\definecolor{wildcolor}{gray}{0.8}


\usepackage{relsize}


%keyword format
\newcommand{\keyw}[1]{{\color{keywordcolor}\texttt{\textbf{#1}}}}

\newcommand{\true}{\emph{true}}
\newcommand{\false}{\emph{false}}


\newcommand\lin{\lstinline[mathescape,basicstyle=\ttfamily]}
\newcommand{\linn}[1] {\begin{center}\lin{1}\end{center}}
\newcommand{\sem}[1]{\llbracket #1\rrbracket}


\renewcommand{\anddot}{\ldotp}

\newcommand{\van}{\abr{WL}\xspace}

\newcommand{\xmark}{\ding{55}}%

%%% THIS PAPER NOTATION
\newcommand{\smt}{\mathcal{L}_\mathcal{T}}
\newcommand{\lrap}{\mathcal{L}_{\text{Pr}(\mathcal{T})}}
\newcommand{\dist}{D}
\newcommand{\pdf}{p}
\newcommand{\pdfs}{\mathcal{D}}

\newcommand{\apprs}{\mathcal{A}}
\newcommand{\appr}{\mathit{step}}
\newcommand{\sappr}{\appr^\phi}
\newcommand{\cdf}{F}
\newcommand{\pop}{\dist_{\emph{pop}}}
\newcommand{\pr}[1]{\text{Pr}[#1]}
\newcommand{\pru}[1]{\widehat{\text{Pr}}[#1]}
\newcommand{\prl}[1]{\widecheck{\text{Pr}}[#1]}

\newcommand{\smtp}[2]{#1|#2}
\newcommand{\var}{\emph{Var}}
\newcommand{\lab}{\emph{lab}}
\newcommand{\expr}{\emph{exp}}
\newcommand{\assume}{\emph{assume}}
\newcommand{\pc}{\emph{pc}}
\newcommand\cond[1]{\color{MidnightBlue}{ \left\{#1 \right\} } }
\newcommand\bcond[1]{\color{MidnightBlue}{\left\{ #1 \right\}}}



%programs
\newcommand{\prog}{\mathcal{P}}
\newcommand{\pre}{\mathcal{P}_\emph{pre}}
\newcommand{\dec}{\mathcal{P}_\emph{dec}}
\newcommand{\post}{\mathit{post}}

\newcommand{\vars}{\mathcal{X}}
\newcommand{\inputs}{\vec{v}_i}
\newcommand{\outputs}{\vec{v}_o}
\newcommand{\stmt}{s}
\newcommand{\stmts}{S}
\newcommand{\bexp}{B}
\newcommand{\aexp}{E}
\newcommand{\sif}{\texttt{\textbf{if}}\xspace}
\newcommand{\sthen}{\texttt{\textbf{then}}\xspace}
\newcommand{\selse}{\texttt{\textbf{else}}\xspace}
\newcommand{\gauss}{\texttt{gauss}\xspace}
\newcommand{\skips}{\ensuremath{\varepsilon}}

\newcommand{\val}{\emph{s}}


% operational semantics
\newcommand{\assign}{\abr{ASSIGN}\xspace}
\newcommand{\passign}{\abr{PASSIGN}\xspace}
\newcommand{\tcond}{\abr{TCOND}\xspace}
\newcommand{\fcond}{\abr{FCOND}\xspace}
\newcommand{\seq}{\abr{SEQ}\xspace}

% operational semantics
\newcommand{\eassign}{\abr{VC-ASN}\xspace}
\newcommand{\epassign}{\abr{VC-PASN}\xspace}
\newcommand{\econd}{\abr{VC-COND}\xspace}
\newcommand{\eseq}{\abr{VC-SEQ}\xspace}


\newcommand{\fs}{FairSquare\xspace}
\newcommand{\wvc}{\abr{wvc}\xspace}
\newcommand{\vol}{\textsc{vol}\xspace}
\newcommand{\volalg}{\abr{SYMVOL}\xspace}
\newcommand{\volopt}{\abr{ADF-SYMVOL}\xspace}
\newcommand{\volume}{\mathit{vol}}
\newcommand{\pair}[1]{\langle #1 \rangle}

\newcommand{\pexp}{\emph{PExp}}
\newcommand{\prob}{\emph{Prob}}

\newcommand{\verify}{\abr{VERIFY}\xspace}

\newcommand{\semprog}{\varphi_\prog}
\newcommand{\enc}{\abr{PVC}\xspace}

\newcommand{\rect}{\mbox{\mancube}}

\newcommand{\hdecomp}{\abr{HDECOMP}\xspace}
\newcommand{\sample}{\abr{HSAMPLE}\xspace}
\newcommand{\decay}{\abr{DECAY}\xspace}

\newcommand{\block}{\mathit{block}}
\newcommand{\rects}{\emph{rects}}

\pgfmathdeclarefunction{gauss}{2}{%
  \pgfmathparse{1/(#2*sqrt(2*pi))*exp(-((x-#1)^2)/(2*#2^2))}%
}


\let\originalleft\left
\let\originalright\right
\renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft}
\renewcommand{\right}{\aftergroup\egroup\originalright}

prelims.tex
\section{Probabilistic Programs and Verification}\label{sec:prelims}

% \begin{table}
%   \centering
% \begin{tabular}{lll}
%   Distribution & \abr{PDF} & \abr{CDF} \\
%   \hline
%   \texttt{gaussian}$(\mu,\sigma)$ &
%   \(\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}\)&
%   0.5\left(1+\text{erf}(\frac{x-\mu}{\sigma \sqrt{2}})\right)
%   \\
%
%   \texttt{uniform}$(a,b)$ &
%   \(\begin{cases} (b-a)^{-1}& x \in[a,b]\\0 & \text{o/w}\end{cases}\) &
%   \(\begin{cases} 0 & x < a\\ \frac{x - a}{b-x} & x\in [a,b) & 1 & x \geq b\end{cases}\) &
%
%    \\
%
%   \texttt{exponential}$(\lambda)$ & \(\lambda e^{-\lambda x}\) & \(1-e^{-\lambda x}\) \\
%   \texttt{laplace}$(\mu, b)$ &
%   \(\(\frac{1}{2b}\text{exp}\left(-\frac{|x-\mu|}{b}\right)\) & \\
%   \ldots&&\\
%   \hline
% \end{tabular}
% \caption{Example univariate probability distributions}
% \end{table}
We formally define programs
and present a general framework for stating and verifying probabilistic
properties.
%
% \begin{figure*}[t]
%   \centering
%   \smaller
%   \begin{prooftree}
%   \val' = \val[x \mapsto \val(e)]
%   \justifies
%   (\omega, \val, x \gets e) \rightarrow (\omega, \val', \skips)
%   \using \assign
%   \end{prooftree}
%   %
%   \hspace{.5in}
%   %
%   \begin{prooftree}
%   \val' = \val[x \mapsto c]
%   \justifies
%   (c:\omega, \val, x \sim p) \rightarrow (\omega, \val', \skips)
%   \using \passign
%   \end{prooftree}
%   %
%   \hspace{.5in}
%   %
%   \begin{prooftree}
%   (\omega, \val, S_1) \rightarrow (\omega', \val', S_1')
%   \justifies
%   (\omega, \val, S_1S_2) \rightarrow (\omega', \val', S_1'S_2)
%   \using \seq
%   \end{prooftree}
%   ~\\
%   \vspace{.1in}
%
%   \begin{prooftree}
%     ~\vspace{.085in}
%   \justifies
%   (\omega, \val, \skips S) \rightarrow (\omega, \val, S)
%   \using \textsc{skip}
%   \end{prooftree}
% %
%   \hspace{.4in}
% %
%   \begin{prooftree}
%   \val(b) = \true
%   \justifies
%   (\omega, \val, \sif~b~\sthen~S_1~\selse~S_2) \rightarrow (\omega, \val, S_1)
%   \using \tcond
%   \end{prooftree}
%   %
%   \hspace{.4in}
%   %
%   \begin{prooftree}
%   \val(b) = \false
%   \justifies
%   (\omega, \val, \sif~b~\sthen~S_1~\selse~S_2) \rightarrow (\omega, \val, S_2)
%   \using \fcond
%   \end{prooftree}
%   \caption{Operational semantics}
%   \label{fig:sem}
% \end{figure*}

\subsection{Program model and  semantics}
\paragraph{Programs}
A program $\prog$ is a sequence of statements $S$:
\begin{align*}
    S \coloneqq & ~V \gets \aexp & \text{ assignment statement}\\[-3pt]
        \mid&~ V \sim \dist &\text{ probabilistic assignment}\\[-3pt]
        \mid&~ \sif~ B ~\sthen~ S ~\selse~ S & \text{conditional}\\[-3pt]
        \mid&~ SS & \text{ sequence of statements}
\end{align*}
%
where $V$ is the set of \emph{real-valued variables} that can appear
in $\prog$, $e \in \aexp$ is an arithmetic expression over variables in $V$,
and $b \in B$ is a Boolean expression over variables in $V$.
A probabilistic assignment is made by sampling from
a probability distribution $p \in D$. A probability distribution
can be, for example,
a \emph{Gaussian distribution}, denoted by $\gauss(\mu, \sigma)$, where
$\mu, \sigma \in \mathds{R}$ are the mean and standard deviation
of the Gaussian.
We shall  restrict distributions to be \emph{univariate} and with
constant parameters, e.g., mean and standard deviation of a Laplacian
or Gaussian---that is, we assume independence of probabilistic assignments.
% \footnote{
% Our language is also able to capture probabilistic assignments of the form
% $y \sim \gauss(x, \sigma)$ where $x\in V$ and $\sigma\in\mathds{R}$.
% This can be transformed into an equivalent
% sequence of assignments $y\sim \gauss(0, \sigma); y\gets y+x$.
% }
Given a probabilistic assignment $x \sim p$, we shall
treat $p(x)$ as a \emph{probability density function} (\abr{PDF})
of the distribution from which the value assigned to $x$ is drawn.
For instance, if the distribution is \texttt{gauss(0,1)},
then $p(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{\emph{x}^2}{2}}$.
%The \emph{cumulative distribution function} (\abr{CDF}) of $p(x)$
%is $F(x) = \int_{-\infty}^x p(x)~ dx$.
%

We use $\vec{v}_i$ to denote a vector of \emph{input variables}
of $\prog$, and $\vec{v}_o$ to denote a vector of \emph{output variables}
of $\prog$; these variables appear in $V$ and
denote the arguments and returns of $\prog$.
We say that a program is \emph{closed} if it has no inputs, i.e.,
$\vec{v}_i$ is empty.
We shall  refer to the following subsets of $V$.
\begin{itemize}
    \item $V_p \subseteq V$ is the set of \emph{probabilistic variables}:
    those that get assigned to in probabilistic assignments.
    \item $V_d = V \setminus V_p$ is the set of \emph{deterministic variables}:
    those that do not appear in probabilistic assignments.
\end{itemize}

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).%
\footnote{Unbounded loops
could be handled through iterative unrolling
as in~\cite{sankaranarayanan13}.}
As demonstrated in Sec.~\ref{sec:example},
the same language is used to define probabilistic preconditions programmatically.


\paragraph{Operational semantics}
% We now define the operational semantics
% of our program model.
Typically,
the \emph{state} $\val : V \rightarrow \mathds{R}$  of the program is defined
as a valuation function from variables in $V$ to
values in $\mathds{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~\cite{kozen1981semantics},
we assume a finite sequence of independent
 random variables.
The semantics of an execution is thus defined for a
fixed sequence of values $\omega$ 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
\iftoggle{long}{Appendix~\ref{app:opsem}}{\emph{supplementary materials}}.
% 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 $\prog$ as shown in Figure~\ref{fig:sem}.
% We use $\val(e)$ to denote the value of expression
% $e$ under a given state $\val$, and use $\skips$
% to denote the empty statement sequence.
% We use the substitution notation
% $s[x \mapsto c]$ to denote the state $s$ but with
% $x$ mapped to $c \in \mathds{R}$.
% The interesting rule is \passign, which, given a probabilistic
% assignment $x \sim \pdf$, for some variable $x$ and distribution $\pdf$,
% picks and removes the first element $c$ of the sequence $c:\omega$
% 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 $\prog$, we define $\sem{\prog}_\omega$
% as the final state  $s$ reachable from executing $\prog$
% using the sequence $\omega$, as defined by the
% relation $\rightarrow$ in Figure~\ref{fig:sem}.

%



\subsection{Programs and volume computation}
% We are often interested in computing the probability
% that a closed program $\prog$ terminates
% in a state satisfying some predicate $\varphi$ over $V$.
% In the case where the program manipulates finite domains,
% e.g., Booleans, the problem reduces to \emph{weighted model counting}~\cite{chavira2008probabilistic},
% where one counts all executions that satisfy $\varphi$
% and weighs them by their probability of occurrence.
% In our setting, programs operate  over the dense field of reals,
% and therefore, instead of counting, we need to \emph{integrate} over the
% infinite space of program executions.

Following Chistikov et al.~\cite{chistikov14arxiv}, we reduce the
problem of computing the probability that the program
terminates in a state satisfying $\varphi$
to \emph{weighted volume computation} (\wvc) over formulas
 describing regions in $\mathds{R}^n$.
In what follows, we begin by formalizing the \wvc problem.

\paragraph{Volume of a formula}
We will use $\smt$ to denote  first-order
formulas in some real arithmetic theory $\mathcal{T}$.
Specifically, we consider two decidable theories:
\emph{linear real arithmetic} and the strictly
richer \emph{real closed fields}---Boolean combinations
of polynomial inequalities.
Given a formula $\varphi \in \smt$, a model $m$ of $\varphi$,
denoted by $m \models \varphi$, is a point in $\mathds{R}^n$,
where $n$ is the number of free variables of $\varphi$.
Thus, we view $\varphi$ as a region in $\mathds{R}^n$, i.e.,
$\varphi \subseteq \mathds{R}^n$.
We use $\vars_\varphi = \{x_1,\ldots,x_n\}$
to denote the free variables
of $\varphi$.

The (unweighted) \emph{volume} of a formula $\varphi$ is
$\int_\varphi 1~ d\vars_\varphi$
where $d\vars_\varphi$ is short for $dx_1dx_2\ldots dx_n$.
For example, if $\varphi$ is in $\mathds{R}^2$, then
$\int_\varphi 1~d\vars_\varphi$ is the area of $\varphi$.


\paragraph{Weighted volume of a formula}
We now define the \emph{weighted volume} of a formula.
We assume we are given a pair $(\varphi, \pdfs)$,
where $\varphi \in \smt$ and $\pdfs = \{\pdf_1,\ldots,\pdf_n\}$ is a set
of  probability density functions
such that each variable $x_i \in \vars_\varphi$ is
associated with a density function $\pdf_i(x_i)$
of the probability distribution of its values.
The weighted volume of $\varphi$ with respect to $\pdfs$, denoted by $\vol(\varphi,\pdfs)$,
is defined as follows:
%
$$\textstyle \int_\varphi \prod_{x_i \in \vars_\varphi} \pdf_i(x_i) ~d\vars_\varphi$$
%
%where $d\vars_\varphi$ is short for $dx_1,\ldots,dx_n$.
\begin{example}
  Consider the formula $\varphi \equiv x_1 + x_2 \geq 0$,
  and let $\pdfs = \{p_1,p_2\}$, where $p_1$ and $p_2$ are
  the \abr{PDF} of the
  Gaussian distribution with mean $0$ and standard deviation $1$.
  Then,
  $$\textstyle \vol(\varphi,\pdfs) = \int_{x_1+x_2\geq0} p_1(x_1)p_2(x_2)~dx_1dx_2 = 0.5$$
  Intuitively, if we are to randomly draw two values for $x_1$ and $x_2$
  from the Gaussian distribution,
  we will land in the region $x_1+x_2\geq 0$ with  probability 0.5.
\end{example}

\paragraph{Probabilistic verification conditions}
Recall that our goal is to compute the probability
of some predicate $\varphi$ at the end of a program
execution, denoted $\pr{\varphi}$.
We now show how to encode this problem
as weighted volume computation.
First, we encode program executions
as a formula $\semprog$.
The process is similar to standard
verification condition generation (as used by verification~\cite{barnett2005weakest}
and bounded model checking tools~\cite{clarke2004tool}), with the difference
that probabilistic assignments populate a set $\pdfs$
of probability density functions.
%

\begin{figure}[t]
  \centering
  \smaller
  \begin{prooftree}
  \justifies
  \pair{x = \sem{e}, \emptyset} \rhd x \gets e
  \using \eassign
  \end{prooftree}
%
%  \vspace{.05in}
  \hspace{1em}
  \begin{prooftree}
  \justifies
  \pair{\true, \{\pdf_i\}} \rhd x_i \sim \pdf_i
  \using \epassign
  \end{prooftree}

  \vspace{.1in}
  \begin{prooftree}
  \pair{\varphi_1,\pdfs_1} \rhd S_1
  \quad\quad
  \pair{\varphi_2,\pdfs_2} \rhd S_2
  \justifies
  \pair{\varphi_1 \land \varphi_2, \pdfs_1 \cup \pdfs_2} \rhd S_1S_2
  \using \eseq
  \end{prooftree}

  \vspace{.1in}
  \begin{prooftree}
  \pair{\varphi_1,\pdfs_1} \rhd S_1
  \quad\quad
  \pair{\varphi_2,\pdfs_2} \rhd S_2
  \justifies
  \pair{\emph{ite}(\sem{b},\varphi_1,\varphi_2),\pdfs_1 \cup \pdfs_2} \rhd \sif~b~\sthen~S_1~\selse~S_2
  \using \econd
  \end{prooftree}
  %
  \caption{Probabilistic verification condition generation.
  Above, $\emph{ite}(a,b,c) \triangleq (a \Rightarrow b) \land (\neg a \Rightarrow c)$.}
  \label{fig:enc}
\end{figure}


 Figure~\ref{fig:enc} inductively defines
 the construction of a \emph{probabilistic verification condition} for a program $\prog$, denoted by a function $\enc(\prog)$, which returns a pair $\pair{\semprog, \pdfs}$.
Without loss of generality, to simplify our exposition,
we assume programs
are in \emph{static single assignment} (\abr{SSA}) form~\cite{cytron1991efficiently}. Given a Boolean expression $b$,
the denotation $\sem{b}$ is the same expression
interpreted as an $\smt$ formula.
The same applies to arithmetic expressions $e$.
For example, $\sem{\texttt{x + y > 0}} \triangleq x + y > 0$.
Intuitively, the construction generates a formula $\varphi_\prog$
that encodes program executions, treating probabilistic assignments
as non-deterministic,
and a set $\pdfs$ of the \abrs{PDF} of distributions in probabilistic
assignments (rule \epassign).

Now, suppose we are given a closed program $\prog$ and a
Boolean formula $\varphi$ over its output variables. Then,
$$\pr{\varphi} = \vol(\exists V_d \ldotp \semprog \land \varphi,\pdfs)$$
That is, we project out all non-probabilistic variables
from $\semprog \land \varphi$ and compute the weighted
volume with respect to the densities $\pdf_i \in \pdfs$.
Intuitively, each model $m$ of $\exists V_d \ldotp \semprog \land \varphi$
corresponds to a sequence of values drawn
in probabilistic assignments in an execution of $\prog$.
%sequence $\omega$ such that $\sem{\prog}_\omega$ satisfies
%$\varphi$.
We note that our construction is closely
related to that of Chistikov et al.~\cite{chistikov14arxiv},
to which we refer the reader for a measure-theoretic
formalization.

\begin{example}
Consider the following closed program $\prog$
\begin{lstlisting}
  x ~ gauss(0,2); y ~ gauss(-1,1); z $\gets$ x + y
\end{lstlisting}
where \texttt{z} is the return variable.
Using the encoding in Figure~\ref{fig:enc},
we compute the pair $\pair{\semprog,\pdfs} \rhd \prog$, where
$\semprog \triangleq z = x + y$
and
$\pdfs = \{\pdf_x, \pdf_y\}$,
where $\pdf_x$ and $\pdf_y$ are the \abrs{PDF} of the
two distributions from which values of $x$ and $y$
are drawn.
% For example, the function
% $\pdf_x(x) = \frac{1}{2\sqrt{2\pi}}e^{-\frac{x^2}{8}$.

Suppose that we would like to compute the probability
that \texttt{z} is positive when the program terminates:
$\pr{z \geq 0}$.
Then, we can compute the following weighted volume:
$\vol(\exists z \ldotp \semprog \land z \geq 0, \pdfs)$,
which is equal to $0.32736$.
%$\int_{x + y \geq 0} \pdf_x(x)\pdf_y(y) ~ dx dy = 0.32736,$
%where $x + y \geq 0$ and
% $\exists z \ldotp z = x + y \land z \geq 0$ are equivalent.
\end{example}

\subsection{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.

\paragraph{Verification problems}
A verification problem is a triple
$(\pre, \dec, \post)$,
where
\begin{itemize}
  \item  $\pre$,  called the \emph{probabilistic precondition},
   is a closed program
  over variables $V^\emph{pre}$ and
  output variables $\outputs^\emph{pre}$.

  \item $\dec$, called the \emph{decision-making program},
  is an open program over variables $V^\emph{dec}$;
  its input arguments are $\inputs^\emph{dec}$,
  with $|\inputs^\emph{dec}| = |\outputs^\emph{pre}|$;
  and its output variables are $\outputs^\emph{dec}$.
  (We assume that $V^\emph{pre}\cap V^\emph{dec}=\emptyset$.)

  \item $\post$ is a \emph{probabilistic postcondition},
  which is a Boolean expression over probabilities of program outcomes.
  Specifically, $\post$ is defined as follows:
  \newcommand{\oper}{X}
  \begin{align*}
    \post \in \pexp \coloneqq&~ \prob > \prob \mid \pexp \lor \pexp\\
          \mid&~ \pexp \land \pexp \mid \neg \pexp\\
    \prob \coloneqq&~ \pr{\varphi} \mid \pr{\varphi \mid \varphi} \mid c
          \mid \prob \odot \prob\\
    \odot \in&~ \{+,-,\div,\times\} ~~~~ c \in \mathds{R}
  \end{align*}
  where $\varphi \in \smt$ are formulas over input and output
  variables of $\dec$.
  For example, \emph{post} might be of the form
  $$\pr{x > 0} > 0.5 \land \pr{y + z > 7} > \pr{t > 5}$$
\end{itemize}
The goal of verification is to prove that $\post$ is true
for the  program $\dec \circ \pre$, i.e.,
the composition of the two programs where we first
run $\pre$ to generate an input for $\dec$.
Since $\pre$ is closed, the program $\dec \circ \pre$ is also closed.

% \paragraph{Fairness verification problems}
% An $\epsilon$-\emph{fairness} verification problem $(\pre, \dec, \post)$ is the
% special class of verification problems
% where $\post$ is of the form:
% $$\frac{\pr{\varphi_f \mid \varphi_s}}{
%         \pr{\varphi_f \mid \neg \varphi_s}} > 1 - \epsilon$$
% where $\epsilon \in [0,1]$;
% formula $\varphi_f$, called the \emph{fairness target},
% is  over $\outputs^\emph{dec}$;
% and formula $\varphi_s$, called the \emph{sensitive condition},
% is  over $\outputs^\emph{pre}$.
% A fairness verification problem is \emph{fair} if and only if the postconditon holds
%  (as defined below).

\paragraph{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~\ref{alg:verify},
takes a verification problem and returns
whether the probabilistic postcondition holds.

\verify begins by encoding the composition of the two programs, $\dec \circ \pre$, as the pair $\pair{\varphi_\prog, \pdfs}$
and adds the constraint
$\inputs^\emph{dec} = \outputs^\emph{pre}$
to connect the outputs of $\pre$ to the inputs of $\dec$
(recall the example from Sec.~\ref{sec:example} for an illustration).
For each term of the form $\pr{\varphi}$ 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{\varphi}$ with their values in $m$---then the
postcondition holds.
% Notice that in this algorithm we are only interested in computing probabilities of program outcomes.
% This aspect allows us to compute the probabilities of individual events appearing in $post$ independently.

\begin{figure}[!t]
  \small
  \begin{spacing}{1}
    \begin{algorithmic}[1]
        \Function{\verify}{$\pre, \dec, \post$}
        \State $\pair{\varphi_\emph{pre},\pdfs_\emph{pre}} \gets \enc(\pre)$
        \State $\pair{\varphi_\emph{dec},\pdfs_\emph{dec}} \gets \enc(\dec)$
        \State $\pair{\varphi_\prog,\pdfs} \gets
          \pair{\varphi_\emph{pre}\land \varphi_\emph{dec} \land \inputs^\emph{dec} = \outputs^\emph{pre}, \pdfs_\emph{pre} \cup \pdfs_\emph{dec}}$
        \State $V_d \gets V_d^\emph{pre} \cup V_d^\emph{dec}$
        \State $m \gets \emptyset$
        \For{$\pr{\varphi} \in \post$}
          \State $m \gets m[\pr{\varphi} \mapsto \vol(\exists V_d \ldotp \semprog \land \varphi, \pdfs)]$
        \EndFor
        \vspace{-3pt} %for some reason the next line has a small offset
        \State \Return $m \models \post$
        \EndFunction
    \end{algorithmic}
    \caption{Abstract verification algorithm}\label{alg:verify}
  \end{spacing}
\end{figure}

%!TEX root = paper.tex

proofs.tex
\section{Operational Semantics}\label{app:opsem}
\paragraph{Operational semantics}
We now define the operational semantics
of our program model.
Typically,
the \emph{state} $\val : V \rightarrow \mathds{R}$  of the program is defined
as a valuation function from variables in $V$ to
values in $\mathds{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~\cite{kozen1981semantics},
we assume a finite sequence of independent and identically
distributed (\emph{iid}) random variables.
The semantics of an execution is thus defined for a
fixed sequence of values $\omega$ 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 $\prog$ as shown in Figure~\ref{fig:sem}.
We use $\val(e)$ to denote the value of expression
$e$ under a given state $\val$, and use $\skips$
to denote the empty statement sequence.
We use the substitution notation
$s[x \mapsto c]$ to denote the state $s$ but with
$x$ mapped to $c \in \mathds{R}$.
The interesting rule is \passign, which, given a probabilistic
assignment $x \sim \pdf$, for some variable $x$ and distribution $\pdf$,
picks and removes the first element $c$ of the sequence $c:\omega$
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 $\prog$, we define $\sem{\prog}_\omega$
as the final state  $s$ reachable from executing $\prog$
using the sequence $\omega$, as defined by the
relation $\rightarrow$ in Figure~\ref{fig:sem}.


\begin{figure*}[t]
  \centering
  \smaller
  \begin{prooftree}
  \val' = \val[x \mapsto \val(e)]
  \justifies
  (\omega, \val, x \gets e) \rightarrow (\omega, \val', \skips)
  \using \assign
  \end{prooftree}
  %
  \hspace{.5in}
  %
  \begin{prooftree}
  \val' = \val[x \mapsto c]
  \justifies
  (c:\omega, \val, x \sim p) \rightarrow (\omega, \val', \skips)
  \using \passign
  \end{prooftree}
  %
  \hspace{.5in}
  %
  \begin{prooftree}
  (\omega, \val, S_1) \rightarrow (\omega', \val', S_1')
  \justifies
  (\omega, \val, S_1S_2) \rightarrow (\omega', \val', S_1'S_2)
  \using \seq
  \end{prooftree}
  ~\\
  \vspace{.1in}

  \begin{prooftree}
    ~\vspace{.085in}
  \justifies
  (\omega, \val, \skips S) \rightarrow (\omega, \val, S)
  \using \textsc{skip}
  \end{prooftree}
%
  \hspace{.4in}
%
  \begin{prooftree}
  \val(b) = \true
  \justifies
  (\omega, \val, \sif~b~\sthen~S_1~\selse~S_2) \rightarrow (\omega, \val, S_1)
  \using \tcond
  \end{prooftree}
  %
  \hspace{.4in}
  %
  \begin{prooftree}
  \val(b) = \false
  \justifies
  (\omega, \val, \sif~b~\sthen~S_1~\selse~S_2) \rightarrow (\omega, \val, S_2)
  \using \fcond
  \end{prooftree}
  \caption{Operational semantics}
  \label{fig:sem}
\end{figure*}

\section{Proofs} \label{app: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 $\rect_\varphi$ characterize all hyperrectangles in $\varphi$
and no others.

\begin{theorem}[Correctness of $\rect$]\label{thm:decomp}
Let $\varphi \in \smt$.
%\begin{itemize}
%\item
  \textbf{Soundness}:
Let $m \models \rect_\varphi$. Then, $H^m \Rightarrow \varphi$
is valid.
%\item
\textbf{Completeness}:
Let $H$ be a hyperrectangle such that  $H \Rightarrow \varphi$.
Then, the following is satisfiable:
$\rect_\varphi \land \bigwedge_{x \in \vars_\varphi} l_x = H_l(x) \land u_x = H_u(x) $
%\end{itemize}
\end{theorem}

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

\begin{theorem}[Correctness of $\block$] \label{thm:block}
Given $\varphi$, let $\Psi \Rightarrow \rect_\varphi$, and
let $m_1,m_2 \models \Psi$.
%\begin{itemize}
%
%\item
\textbf{Soundness}:
    If $H^{m_1} \land H^{m_2}$ is satisfiable,
    then $m_2 \not\models \Psi \land \block(H^{m_1})$.
%
%\item
    \textbf{Completeness}: If $H^{m_1} \land H^{m_2}$
is unsatisfiable,
then $m_2 \models \Psi \land \block(H^{m_1})$.
%\end{itemize}
\end{theorem}

The following theorem states the correctness of the \abr{ADF}
encoding:

\begin{theorem}[Correctness of $\sappr$]\label{thm:adf}
Fix an \abr{ADF} $\appr(x)$.
%\begin{itemize}
  %\item
  \textbf{Soundness}:
For any model $m \models \sappr$,
the following  is true:
$m(\delta_x) = \int_{m(l_x)}^{m(u_x)} \appr(x)~dx$.
%\item
\textbf{Completeness}:
For any constants $a,b,c\in \mathds{R}$ such that
$c = \int_a^b \appr(x)~dx$, the following formula is satisfiable:
$\delta_x = c \land l_x = a \land u_x = b \land \sappr(x)$.
%\end{itemize}
\end{theorem}


\subsection*{Proof of Theorem~\ref{thm:decomp}}
\paragraph{Soundness} Suppose $m\models \rect_\varphi$.
By definition of $\rect_\varphi$,
the following formula is valid $$\bigwedge_{x \in \vars_\varphi} m(l_x) \leq x \leq m(u_x) \Rightarrow \varphi$$
Therefore $H^m \Rightarrow \varphi$, by definition of $H^m$.

\paragraph{Completeness}
Let $H$ be a hyperrectangle such that $H \Rightarrow \varphi$.
By definition,
$H$ is of the form $\bigwedge_{x\in \vars_\varphi} c_x \leq x \leq c_x'$.
It immediately follows that the model where $l_x = c_x$ and $u_x = c_x'$,
for every $x\in \vars_\varphi$, satisfies $\rect_\varphi$,
since $c_x' \geq c_x$ (satisfying the first conjunct of $\rect_\varphi$),
and $\forall \vars_\varphi \ldotp H \Rightarrow \varphi$ (satisfying the second conjunct of $\rect_\varphi$).

\subsection*{Proof of Theorem~\ref{thm:block}}
\paragraph{Soundness}
Suppose $H^{m_1} \land H^{m_2}$ is satisfiable.
By definition of a hyperrectangle,
this means that for all  variables $x\in \vars_\varphi$,
we have that the intervals $[H_l^{m_1}(x), H_u^{m_1}(x)]$
and $[H_l^{m_2}(x), H_u^{m_2}(x)]$ overlap, i.e., at least are equal
on one of the extremes.
Therefore, $m_2 \not\models \Psi \land \block(H^{m_1})$,
since $\block$ does not admit any model $m$ where,
for all $x\in\vars_\varphi$,
$[m(l_x),m(u_x)]$ overlaps with $[H_l^{m_1}(x), H_u^{m_1}(x)]$.

\paragraph{Completeness}
Suppose $H^{m_1} \land H^{m_2}$ is unsatisfiable.
By definition of a hyperrectangle,
there is at least one $x \in \vars_\varphi$
where $[H_l^{m_1}(x), H_u^{m_1}(x)]$
and $[H_l^{m_2}(x), H_u^{m_2}(x)]$ do not overlap.
Therefore, if $m_2 \models \Psi$, then $m_2 \models \Psi \land \block(H^{m_1})$,
since $\block$ explicitly states that for at least one variable $x$, $[m(l_x),m(u_x)]$ should not overlap with $[H_l^{m_1}(x), H_u^{m_1}(x)]$.

\subsection*{Proof of Theorem~\ref{thm:volsound}}
At any point in the execution,
$\volume = \sum_{i=1}^l \int_{H_i} \prod p_i(x_i) ~dx_i$,
where $l$ is the number of applications of $\sample$
and $H_i$ is the hyperrectangle sampled at step $i$.
By definition, $\bigvee{H_i} \Rightarrow \varphi$.
Since \abrs{PDF} are positive functions,
$\volume \leq \vol(\varphi,\pdfs)$.

\subsection*{Proof of Corollary~\ref{col:upper}}
By definition of \abrs{PDF} and integration,
$$\int_{\mathds{R}^n} \prod p_i(x_i) ~dx_i =
\int_\varphi \prod p_i(x_i) ~dx_i + \int_{\neg\varphi} \prod p_i(x_i) ~dx_i$$
for any $\varphi \subseteq \mathds{R}^n$.
From Theorem~\ref{thm:volsound}, it follows that
at any point in the execution of $\volalg(\neg\varphi,\pdfs)$,
we have $1 - \volume \geq \vol(\varphi,\pdfs)$.

\subsection*{Proof of Theorem~\ref{thm:adf}}
\paragraph{Soundness}
Suppose $m\models \sappr(x)$.
Then,
$$m(\delta_x) =  \sum_{i = 1}^n c_i \cdot \bigl|[a_i,b_i] \cap [m(l_x),m(u_x)]\bigr|$$
%It is easy to see that $\sqcap$ computes the length of the two intervals.
By definition of the area under a positive step function,
we have $$m(\delta_x) = \int_{m(l_x)}^{m(u_x)} \appr(x)~dx$$

\paragraph{Completeness}
Completeness easily follows from correctness of the encoding
of integrals over step functions as sums.

\subsection*{Proof of Theorem~\ref{thm:limit}}
The algorithm constructs two series in parallel:
the actual volume computation series $\sum v_i$
and the approximated series $\sum a_i$,
where each $v_i$ and $a_i$ 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
$$v_i^\Sigma = \sum\limits_{j=1}^i v_j ~~~~~~~~ a_i^\Sigma = \sum\limits_{j=1}^i a_j$$
It is maintained that
\begin{align*}
    \forall i\ldotp v_i^\Sigma \leq \emph{EVol}_{R \cap \varphi} =& \int_{R \cap \varphi} \prod \pdf(x)~ d\vars_\varphi\\
    \forall i\ldotp a_i^\Sigma \leq \emph{AVol} =& \int_{R \cap \varphi} \prod \appr(x) ~d\vars_\varphi
\end{align*}

Suppose, for the sake of obtaining a contradiction, that our sequence of samples to construct
$\{v_i^\Sigma\}$ and $\{a_i^\Sigma\}$
never, in the limit, samples some subregion $R' \subsetneq R$,
but $H \subseteq R \setminus R'$ is a hypercube such that
$\int_H \prod\pdf(x) ~d\vars_\varphi$ is non-zero.
Then, $\{v_i^\Sigma\}$ approaches some limit that is at most
$v^\Sigma \leq \emph{EVol}_{R \cap \varphi} - \int_H \prod \pdf(x) ~d\vars_\varphi$.
Let $\delta = \int_H \prod \appr(x) ~d\vars_\varphi$.
The sequence $\{a_i^\Sigma\}$ is  monotonically increasing and bounded above by $\emph{AVol}$:
it follows that $\{a_i^\Sigma\}$ converges to some limit $a^\Sigma \leq \emph{AVol}$.
By the definition of a limit, for all $\epsilon > 0$, there exists $N$ such that
for all $n > N$, $a^\Sigma - a_n^\Sigma < \epsilon$.
So at some point when we have fixed a threshold $\tau < \delta$
and have run out of samples in $R'$ with
$a_n \geq \tau$ (guaranteed when $a^\Sigma - a_n^\Sigma < \tau$)
we would sample $H \subseteq R \setminus R'$.
This property ensures that in the limit, $R' \rightarrow R$, and
$v_i^\Sigma \rightarrow \emph{EVol}_{R \cap \varphi}$.

\subsection*{Proof of Theorem~\ref{thm:limitfull}}
More formally, $\volopt$ using each $P_i$ and its $\apprs_i$
infinitely often in the fair serialization
translates to the existence of a surjective function $f:\mathds{N}\rightarrow\mathds{N}$
such that \rone the $n$th successful call to $\sample$ is in $P_{f(n)}$
and \rtwo for all $n'$ the preimage $f^{-1}(n')$ is an infinite set.

Since the domain ($R_i$) of each $\apprs_i$ is disjoint from the others,
since $\bigcup_{i=1}^\infty R_i = \mathds{R}^n$,
and since each $P_i$ converges to the weighted volume over its domain,
it follows that
$$\sum_{i=1}^\infty \sum_{j=1}^\infty \vol(H_{i j},\pdfs) = \vol(\varphi,\pdfs)$$
where $H_{i j}$ is the $j$th hyperrectangle returned from $P_i$.
(The $n$th hyperrectangle produced from the serialization corresponds to
$H_{i j}$ with $i = f(n)$ and $j = |f^{-1}(i) \cap \{n' \mid n' \leq n\}|$.)

In fact, since all the terms are non-negative,
the series above \emph{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
$$\sum_{k=1}^\infty \vol(H_k,\pdfs) = \vol(\varphi,\pdfs)$$
Thus, by composing rearrangements, any fair serialization also converges.

% \subsection*{Proof of Theorem~\ref{thm:fair}}
% We assume that division by zero is $\infty$
% and that $\infty > c$, for all $c \in \mathds{R}$.
% \paragraph{Soundness}
% It follows from soundness of the verification condition generation
% and $\volalg$ that $$\pru{\varphi} \geq \pr{\varphi} \geq \prl{\varphi}$$
% for any program property $\varphi$.
% Assume the algorithm returns \emph{fair}. Then,  we know that
% $$\frac{\prl{\varphi_f \land \varphi_s} * \prl{\neg \varphi_s}}
% {\pru{\varphi_f \land \neg \varphi_s} * \pru{\varphi_s}}
% > 1-\epsilon$$
% By definition of upper/lower bounds, we know that
% $$\frac{\prl{\varphi_f \land \varphi_s} * \prl{\neg \varphi_s}}
% {\pru{\varphi_f \land \neg \varphi_s} * \pru{\varphi_s}}
% \leq
% \frac{\pr{\varphi_f \land \varphi_s} * \pr{\neg \varphi_s}}
% {\pr{\varphi_f \land \neg \varphi_s} * \pr{\varphi_s}}$$
% Therefore, the problem is fair.
%
% An analogous argument applies for the \emph{unfair} case.
%
% \paragraph{Almost completeness}
% Suppose the problem is fair and that the exact value of the
%  ratio is $r \in \mathds{R}$.
%  Then, by our convergence assumption on volume computation,
%  we will eventually reach a lower bound $r'$, where  $1-\epsilon < r' < r$,
%  thus proving fairness.
%
% Suppose the problem is unfair and the exact value of the ratio
% is $r \neq 1-\epsilon$.
% Then, by our convergence assumption, we will eventually
% reach an upper bound $r'$, where $1-\epsilon > r' > r$, thus proving unfairness.

related.tex
\section{Related Work}\label{sec:relwork}
%
% Our work touches upon a very broad spectrum of topics.
% Here, we compare with the most related works.


\paragraph{Probabilistic program analysis}
% There is a rich and diverse body
% of work on analysis of probabilistic programs.
We refer the reader to Gordon et al. for a thorough survey~\cite{gordon2014probabilistic}.
% Probabilistic program semantics were first explored by Kozen~\cite{kozen1981semantics} and later expanded upon~\cite{cousot2012probabilistic,monniaux2005abstract}.
%
% ABSTRACT INTERP
A number of works tackled analysis of probabilistic programs
from an abstract interpretation perspective~\cite{monniaux2001backwards,monniaux2000abstract,monniaux2001abstract,mardziel2011dynamic,claret2013bayesian}.
The comparison between our solution through volume
computation and abstract interpretation is perhaps
analogous to \abr{SMT} solving and software model checking
versus abstract interpretation.
For example, techniques like Monniaux's~\cite{monniaux2000abstract}, 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.

% McIver and Morgan~\cite{mciver2006abstraction}
% pioneered \emph{quantitative invariants},
% a line of work that formalized inference
% of \emph{expectations} of quantities at loop boundaries.
% Katoen et al.~\cite{katoen2010linear} and Chakarov and Sankaranaryanan~\cite{chakarov2014expectation,chakarov2013probabilistic}
% proposed techniques for automatic inference of quantitative invariants.
% While we do not address loops in this paper,
% our approach can naturally handle loops with a fixed
% number of iterations through unrolling, as is common in probabilistic models.
% Additionally, we can potentially extend our approach to unbounded
% loops through iterative unrolling~\cite{sankaranarayanan13}.

% Symbolic execution
A number of works have also addressed probabilistic
analysis through symbolic execution~\cite{filieri2013reliability,sankaranarayanan13,geldenhuys2012probabilistic,sampson2014expressing}.
Filieri et al.~\cite{filieri2013reliability}
and Geldenhuys et al.~\cite{geldenhuys2012probabilistic}
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.




% MCMC
A number of works rely on \emph{sampling} to approximate
the probability of a given program property. The Church~\cite{goodman08}
programming language, for instance, employs the
\emph{Metropolis--Hastings} algorithm~\cite{chib1995understanding},
a \emph{Markov Chain Monte Carlo} (\abr{MCMC}) technique.
Other  techniques perform probabilistic inference
by compiling programs or program paths
to Bayesian networks~\cite{koller2009probabilistic}
and applying \emph{belief propagation}~\cite{minka2012infer} or sampling~\cite{sampson2014expressing} on the network.

% model checking
% Numerous works in model checking have addressed probabilistic, finite-state systems.
% Notably, the \abr{PRISM} tool is a \emph{probabilistic model checker}
% that incorporates algorithms for exact analysis of decidable
% models such as Markov chains and
% probabilistic automata~\cite{kwiatkowska2011prism}.
% These models are finite-state machines that can capture limited classes of
% programs; in particular, \abr{PRISM}  does not handle
% arbitrary probability distributions.
% The population models in this paper are outside the scope
% of models supported by \abr{PRISM}.

% computer algebra
% Computer algebra techniques have also
% been used to characterize the joint posterior distribution of a
% probabilistic program~\cite{Gehr16,NarayananCRSZ16}.
% Even if a closed form exists, such
% techniques are not guaranteed to find it,
% and might resort to approximation~\cite{NarayananCRSZ16}, e.g., \abr{MCMC}.
% Our volume computation is based on  numerical integration and
% is \emph{guaranteed} to converge to the exact solution in our program model.

\paragraph{Volume computation}
%
The computation of weighted volume is known to be hard---even
for a polytope, volume computation is \#P-hard~\cite{khachiyan1993complexity}.
Two general approaches exist:
approximate and exact solutions.
Note that in general, any approximate technique
at best can prove facts \emph{with high probability}.
% Since our interest is in computing a proof of fairness,
% which involves computing probabilities of a number of events,
% exact solutions are a natural choice.

Our volume computation algorithm is inspired
by \rone Li et al.'s~\cite{li2014management}
formula decomposition procedure,
where quantifier elimination is  used
to under-approximate an \abr{LRA}
constraint as a Boolean combination of monadic
predicates;
and \rtwo Sankaranarayanan et al.'s~\cite{sankaranarayanan13}
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~\cite{Bournez99}.)

A number of factors differentiate our work from~\cite{sankaranarayanan13},
which we compared with experimentally in Section~\ref{sec:impl}.
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.
%In practice, we attempted to use the
%volume computation tool of Sankaranarayanan et al. (as acquired from
%the authors);
%unfortunately, even for volume computation problems
%with a single, simple polyhedron, the tool provides
%impractical probability bounds.
%For example, on \svm{4} with the independent population model, an execution of
%their tool took 4 minutes to produce bounds of $[0.039, 0.233]$
%for the quantity $\pr{\varphi_f \land \varphi_s}$;
%longer executions did not appreciably improve these bounds.
%In contrast, \fs reaches a bound of $[0.113,0.121]$ within seconds.
%We believe this is due to the fact that the \emph{ray-shooting}
%approach it employs does not attempt
%to find the  hyperrectangles with largest volume, like we do here,
%and thus ends up sampling hyperrectangles with negligible weighted volume.
%Since the main goal of this paper is proving algorithmic fairness,
%we do not quantitatively compare our volume computation algorithm
%against other tools.

\emph{LattE} is a tool that performs exact integration
of polynomial functions over polytopes~\cite{de2012software}.
Belle et al.~\cite{belle15,belle16}
compute the volume of a linear real arithmetic (\abr{LRA}) formula
by, effectively, decomposing it into \abr{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 \rone  handles formulas over real closed fields, which subsumes \abr{LRA},
and \rtwo handles arbitrary probability distributions.


Chistikov et al.~\cite{chistikov14arxiv,chistikov15}
present a framework for approximate counting with probabilistic
guarantees in \abr{SMT} theories,
which they specialize for bounded \abr{LRA}.
In contrast, our technique
\rone handles unbounded formulas in \abr{LRA}
as well as real closed fields,
\rtwo handles arbitrary distributions,
and \rthree 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~\cite{dyer1991random}
(see Vempala~\cite{vempala2005geometric}
for a survey).

\paragraph{Algorithmic fairness}
Our work is inspired by recent concern in the
fairness of modern decision-making programs~\cite{zarsky2014understanding,barocas2014big}.
A number of recent works have explored algorithmic fairness~\cite{zemel13,feldman15,hardt16,dwork12,calders10,pedreshi2008discrimination,datta2016algorithmic,datta2015automated}.
%
For instance, Zemel et al.~\cite{zemel13} and Feldman et al.~\cite{feldman15}
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).
% Dwork, et al.~\cite{dwork12} investigated \emph{individual fairness}: that
% similar individuals are treated similarly.
% Other works have explored discrimination from a data mining and
%  transformation perspective~\cite{pedreshi2008discrimination}.
%
%The choice of metric affects the outcome of fairness,
%so in this sense the metric is similar to our notion of a population model.

Discrimination in \emph{black-box} systems
has been studied through the lens of statistical analysis~\cite{sweeney2013discrimination,datta2015automated,datta2016algorithmic}.
Notably, Datta et al.~\cite{datta2015automated}
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.

\paragraph{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.,~\cite{dwork12});
further refinements to the use of \abrs{ADF},
e.g. dynamically adapting them during \volopt;
and improving the weighted volume computation algorithm
for better scalability.




% Sankaranarayanan, et al. compute bounds on
% the weighted volume of polytopes through
% hyperrectangle approximations as well~\cite{sankaranarayanan13}.
% Overapproximations are found by dividing each polytope
% into sections and using linear programming solvers
% to find a bounding box for each section.
% An underapproximating hyperrectangle is found differently:
% a point is sampled inside the polytope and expanded along
% each dimension to form an inscribed hyperrectangle.
% While this latter process is quite similar to our sample maximization technique,
% we were able to produce rudimentary examples
% in which it was insufficient for our desired precision.


%!TEX root = paper.tex

results.tex
\section{Quantitative results}\label{app:res}

The full quantitative results of our experiments (Section~\ref{sec:impl})
are shown in table~\ref{tbl:full};
note that the upper-left quadrant of table~\ref{tbl:full}
is simply the left table from figure~\ref{fig:summary}.
The additional rows are the instances
of qualified group fairness where qualification
is true if and only if $\textit{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 \fs proved to be qualitatively very similar.

The second set of columns detail the results
of running Sankaranarayanan et al.'s tool~\cite{sankaranarayanan13}, \abr{VC}.
For the \emph{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 \emph{Res}.
\emph{T} denotes the amount of time (s) the tool ran for;
\abr{TO} denotes that it timed out after 900s without producing a result.

The third set of columns detail the results
of running \abr{PSI} on the benchmarks.
$T$ denotes the amount of time (s) the tool ran for;
again, \abr{TO} denotes that it timed out after 900s without producing a result.
The \emph{Res} column depicts a check mark if \abr{PSI}
terminated with a closed form \abr{CDF},
and an integral sign if it returned
a function containing unevaluated integrals.

\begin{table*}[t]
\centering
\scriptsize
\setlength{\tabcolsep}{2pt}
\newcommand{\fair}[3]{\multicolumn{1}{c}{\cellcolor{green!20}\ding{51}} & {#1} & {#2} & {#3}}
\newcommand{\unfair}[3]{\multicolumn{1}{c}{\cellcolor{red!25}\ding{55}} & {#1} & {#2} & {#3}}
\newcommand{\stime}[4]{{\tiny $\overset{#1}{#2}$} & {#3} & {\abr{TO}} & {#4}}
\newcommand{\qtime}{\abr{TO}$_q$ & - & - & {\abr{TO}}}
\newcommand{\tblrow}[5]{#1 & #2 & #3 & #4 & #5}
\newcommand{\hhlinesep}{\hhline{|=##=##===##===##===|}}
\newcommand{\trih}{{\it Res} & {\it \#} & {\it Vol} & {\it QE}}
\newcommand{\psih}{{\it Res} & {\it T}}
\newcommand{\vch}{{\it Res} & {\it T}}
\newcommand{\crow}[3]{#1 & #2 & #3}
\newcommand{\vcf}[1]{\multicolumn{1}{c}{\cellcolor{green!20}\ding{51}} & {#1}}
\newcommand{\vcu}[1]{\multicolumn{1}{c}{\cellcolor{red!25}\ding{55}} & {#1}}
\newcommand{\vci}[3]{{\tiny $\overset{#1}{#2}$} & {#3}}
\newcommand{\vct}{- & {\abr{TO}}}
\newcommand{\vcc}{- & {\abr{Cr}}}
\newcommand{\ps}[1]{\ding{52} & {#1}}
\newcommand{\pf}[1]{{$\int$} & {#1}}
\newcommand{\pt}{- & {\abr{TO}}}
\newcommand{\thline}{\specialrule{.15em}{.00em}{.00em}}
\renewcommand{\thline}{\hline}
\newcommand{\qdt}[1]{\abr{DT}$_{#1}$ \abr{Q}}
\newcommand{\qdta}[1]{\abr{DT}$_{#1}^{\alpha}$ \abr{Q}}
\newcommand{\qsvm}[1]{\abr{SVM}$_{#1}$ \abr{Q}}
\newcommand{\qsvma}[1]{\abr{SVM}$_{#1}^{\alpha}$ \abr{Q}}
\newcommand{\qnn}[2]{\abr{NN}$_{#1,#2}$ \abr{Q}}
\newcommand{\fsl}{\cline{1-14}}
{\def\arraystretch{1.0}
\begin{tabular}{*{26}{c}}
    \toprule
    \multirow{4}{1.3cm}{\centering \textbf{Decision program}} & \multirow{4}{*}{\bf Acc} &  \multicolumn{12}{c}{\bf \fs} & \multicolumn{6}{c}{\bf \abr{VC}~\cite{sankaranarayanan13}} & \multicolumn{6}{c}{\bf \abr{PSI}~\cite{Gehr16}}\\
    \cmidrule(lr){3-14}
    \cmidrule(lr){15-20}
    \cmidrule(lr){21-26}
    & & \multicolumn{4}{c}{Independent} & \multicolumn{4}{c}{Bayes Net 1} &\multicolumn{4}{c}{Bayes Net 2} & \multicolumn{2}{c}{Ind} & \multicolumn{2}{c}{BN1} & \multicolumn{2}{c}{BN2} & \multicolumn{2}{c}{Ind} & \multicolumn{2}{c}{BN1} & \multicolumn{2}{c}{BN2}  \\
    \cmidrule(lr){3-6}
    \cmidrule(lr){7-10}
    \cmidrule(lr){11-14}
    \cmidrule(lr){15-16} \cmidrule(lr){17-18} \cmidrule(lr){19-20}
    \cmidrule(lr){21-22} \cmidrule(lr){23-24} \cmidrule(lr){25-26}
    & & \trih & \trih & \trih & \psih & \psih & \psih & \vch & \vch & \vch\\
    \midrule
    \tblrow{\dt{4}}{0.79}{\fair{10}{1.3}{0.5}}{\unfair{12}{2.2}{0.9}}{\unfair{18}{6.6}{2.2}}
        &\crow{\vcf{4.8}}{\vcu{3.7}}{\vci{0.07}{3.82}{16.9}}
        &\crow{\ps{0.07}}{\ps{0.22}}{\pf{1.8}}\\
    \tblrow{\dt{14}}{0.71}{\fair{20}{4.2}{1.4}}{\fair{38}{52.3}{11.4}}{\fair{73}{130.9}{33.6}}
        &\crow{\vci{0.76}{1.34}{12.3}}{\vci{0.74}{1.05}{24.2}}{\vci{0.18}{7.74}{81.2}}
        &\crow{\ps{0.16}}{\ps{6.7}}{\pt}\\
    \tblrow{\dt{16}}{0.79}{\fair{21}{7.7}{2.0}}{\unfair{22}{15.3}{6.3}}{\unfair{22}{38.2}{14.3}} 
        &\crow{\vcf{18.7}}{\vcu{22.3}}{\vci{0.07}{5.06}{78.45}}
        &\crow{\ps{0.2}}{\ps{136.4}}{\pt}\\
    \tblrow{\dta{16}}{0.76}{\fair{18}{5.1}{3.0}}{\fair{34}{32.0}{8.2}}{\fair{40}{91.0}{19.4}}
        &\crow{\vcf{23.2}}{\vcf{36.0}}{\vci{0.13}{6.00}{99.8}}
        &\crow{\ps{0.6}}{\pt}{\pt}\\
    \tblrow{\dt{44}}{0.82}{\fair{55}{63.5}{9.8}}{\unfair{113}{178.9}{94.3}}{\unfair{406}{484.0}{222.4}}
        &\crow{\vcf{40.3}}{\vci{0.56}{1.00}{56.3}}{\vci{0.10}{6.02}{214.6}}
        &\crow{\ps{539.3}}{\pt}{\pt}\\ \thline
    \tblrow{\svm{3}}{0.79}{\fair{10}{2.6}{0.6}}{\unfair{10}{3.7}{1.7}}{\unfair{10}{10.8}{6.2}}
        &\crow{\vct}{\vci{0.00}{153.8}{295.8}}{\vct}
        &\crow{\ps{0.06}}{\pf{0.71}}{\pf{4.2}}\\
    \tblrow{\svm{4}}{0.79}{\fair{10}{2.7}{0.8}}{\unfair{18}{13.3}{3.1}}{\unfair{14}{33.7}{20.1}}
        &\crow{\vct}{\vct}{\vct}
        &\crow{\ps{0.09}}{\pf{1.9}}{\pf{8.2}}\\
    \tblrow{\svma{4}}{0.78}{\fair{10}{3.0}{0.8}}{\fair{22}{15.7}{3.2}}{\fair{14}{33.4}{63.2}}
        &\crow{\vct}{\vct}{\vct}
        &\crow{\ps{0.09}}{\pf{1.9}}{\pf{8.3}}\\
    \tblrow{\svm{5}}{0.79}{\fair{10}{8.5}{1.3}}{\unfair{10}{12.2}{6.3}}{\qtime}
        &\crow{\vct}{\vct}{\vct}
        &\crow{\ps{0.14}}{\pf{5.0}}{\pf{18.2}}\\
    \tblrow{\svm{6}}{0.79}{\stime{0.02}{35.3}{634}{2.4}}{\stime{0.09}{3.03}{434}{12.8}}{\qtime}
        &\crow{\vct}{\vct}{\vci{0.00}{10^{12}}{248.0}}
        &\crow{\ps{0.24}}{\pf{57.8}}{\pf{30.1}}\\ \thline
    \tblrow{\nn{2}{1}}{0.65}{\fair{78}{21.6}{0.8}}{\fair{466}{456.1}{3.4}}{\fair{154}{132.9}{7.2}}
        &\crow{\vci{0.00}{691.3}{27.4}}{\vci{0.00}{410.1}{42.3}}{\vci{0.00}{67.1}{80.7}}
        &\crow{\pf{1.1}}{\pf{5.6}}{\pf{22.4}}\\
    \tblrow{\nn{2}{2}}{0.67}{\fair{62}{27.8}{2.0}}{\fair{238}{236.5}{7.2}}{\fair{174}{233.5}{18.2}}
        &\crow{\vci{0.00}{1365}{97.1}}{\vci{0.00}{162.7}{65.44}}{\vci{0.01}{29.2}{107.1}}
        &\crow{\pf{541.8}}{\pf{712.2}}{\pt}\\
    \tblrow{\nn{3}{2}}{0.74}{\stime{0.03}{674.7}{442}{10.0}}{\stime{0.00}{5.24}{34}{55.9}}{\qtime}
        &\crow{\vct}{\vct}{\vct}
    &\crow{\pf{663.7}}{\pt}{\pt}\\ \thline
    \tblrow{\qdt{4}}{}{\fair{12}{2.0}{1.1}}{\unfair{25}{5.6}{1.8}}{\unfair{23}{10.4}{4.2}} \\
    \tblrow{\qdt{14}}{}{\fair{24}{5.7}{2.4}}{\fair{47}{57.1}{13.8}}{\fair{80}{57.1}{13.8}} \\
    \tblrow{\qdt{16}}{}{\fair{27}{10.1}{3.6}}{\unfair{38}{48.9}{15.3}}{\unfair{128}{179.7}{32.6}} \\
    \tblrow{\qdta{16}}{}{\fair{16}{5.3}{3.7}}{\fair{45}{46.2}{9.3}}{\fair{340}{337.9}{44.1}} \\
    \tblrow{\qdt{44}}{}{\fair{50}{67.5}{13.1}}{\unfair{271}{347.5}{285.9}}{\qtime} \\ \fsl
    \tblrow{\qsvm{3}}{}{\fair{22}{3.1}{3.0}}{\unfair{7.9}{27}{4.8}}{\unfair{24}{18.3}{7.4}} \\
    \tblrow{\qsvm{4}}{}{\fair{18}{4.4}{3.9}}{\unfair{31}{42.9}{19.6}}{\unfair{24}{34.8}{148.3}} \\
    \tblrow{\qsvma{4}}{}{\fair{14}{4.2}{1.9}}{\stime{0.82}{1.12}{705}{5.9}}{\fair{48}{94.9}{49.4}} \\
    \tblrow{\qsvm{5}}{}{\fair{18}{6.4}{5.3}}{\unfair{39}{83.2}{54.0}}{\qtime} \\
    \tblrow{\qsvm{6}}{}{\stime{0.24}{31.5}{686}{3.6}}{\stime{0.04}{3.04}{440}{133.4}}{\qtime} \\ \fsl
    \tblrow{\qnn{2}{1}}{}{\fair{50}{11.5}{1.2}}{\fair{107}{57.5}{3.4}}{\stime{0.71}{1.09}{710}{11.8}} \\
    \tblrow{\qnn{2}{2}}{}{\fair{50}{27.4}{3.1}}{\fair{584}{799.6}{13.6}}{\fair{431}{558.5}{21.5}} \\
    \tblrow{\qnn{3}{2}}{}{\stime{0.25}{47.4}{394}{14.3}}{\stime{0.07}{10^{10}}{93}{53.4}}{\qtime} \\
\end{tabular}}
\caption{Benchmarks}
\label{tbl:full}
\end{table*}

vol.tex
\section{Symbolic Weighted Volume Computation}\label{sec:vol}

In this section, we describe our weighted volume computation
algorithm. Recall that, given a formula  $\varphi$
and a set $\pdfs$ defining the \abrs{PDF}
of the distributions of free variables,
our goal is to evaluate the
integral
%
$\textstyle \int_\varphi \prod_{x_i\in \vars_\varphi}
\pdf_i(x_i) ~d\vars_\varphi$.
%
% \todo{I thought the core challenge was because of the probabilities.}
% The core challenge in solving the above
% problem is the fact the region $\varphi$
% is an arbitrary Boolean formula.
% For example, if $\varphi$ is in linear real arithmetic,
% then it describes potentially exponentially many polyhedra \todo{exponential in what?
% Also, our algorithms doesn't do necessarily better than the number of regions}
% in $\mathds{R}^n$.

\paragraph{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~\cite{dyer1988complexity}.
Existing techniques suffer from
one or more of the following: they
\rone  restrict $\varphi$
to a conjunction of linear inequalities~\cite{sankaranarayanan13,de2012software},
\rtwo restrict integrands to polynomials or simple distributions~\cite{de2012software,belle15,chistikov15,belle2015hashing},
\rthree compute approximate solutions with probabilistic guarantees~\cite{chistikov14arxiv,chistikov15,vempala2005geometric,belle2015hashing},
\rfour restrict $\varphi$ to bounded regions of $\mathds{R}^n$~\cite{chistikov14arxiv,chistikov15},
or \rfive have no convergence guarantees, e.g., computer algebra
tools that find closed-form solutions~\cite{matlab,mathematica,Gehr16}.
(See Sec.~\ref{sec:relwork} for details.)

\paragraph{Symbolic weighted volume computation}
% In this section, we present a general symbolic algorithm
% for solving such integration problems.
Our approach is novel in its generality and its algorithmic core.
The following are the high-level properties of our algorithm:
\begin{enumerate}
  \item Our approach accepts formulas in the decidable yet rich
 theory of \emph{real closed fields}:
  Boolean combinations of polynomial inequalities.

  \item Our approach imposes no restrictions on the form of the
  \abrs{PDF}, only
  that we can evaluate the \emph{cumulative distribution
  functions}\footnote{The \emph{cumulative distribution function} of a real-valued random variable $X$ is the function
  $f:\mathds{R}\rightarrow\mathds{R}$, such that $f(x)=\pr{X\leq x}$.}
  (\abrs{CDF}) associated with the \abrs{PDF}
  in $\pdfs$.

  \item 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.
\end{enumerate}

At the algorithmic level, our approach makes the following contributions:
\begin{enumerate}
  \item Our approach exploits the power of \abr{SMT}
  solvers and uses them as a black box,
  allowing it to directly benefit from future advances in solver technology.
  \item 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~\cite{Bournez99, sankaranarayanan13,asarin2000approximate,li2014management},
  we utilize it in a new symbolic way to enable volume computation over
  \abr{SMT} formulas.
  \item Our approach introduces a novel technique for approximately encoding
  \abrs{PDF} as formulas, and using them to guide the \abr{SMT}
  solver towards making large leaps to the exact solution.
\end{enumerate}

% The only assumptions we make are
% \rone that formulas in the decidable first-order theory of \emph{real closed fields},
% and \rtwo that we can evaluate the \emph{cumulative distribution
% functions}\footnote{The \emph{cumulative distribution function} of a real-valued random variable $X$ is a function
% $f:\mathds{R}\rightarrow\mathds{R}$, such that for any $x\in \mathds{R}$, $f(x)=\pr{X\leq x}$ is the probability that $X$
% takes a value less than or equal to $x$.}
% (\abrs{CDF}) associated with the \abrs{PDF}
% in $\pdfs$.
% Our algorithm monotonically grows the solution
% starting from $0$ and maintains the invariant
% that the current solution is a lower bound of the exact solution.
% We prove that the algorithm converges to the exact solution in the limit.
% An advantage of our technique is that it reduces the problem
% to a series of satisfiability queries to  a
% black-box \abr{SMT} solver, allowing
% us to directly exploit future advances in \abr{SMT} technology.
% %and state conditions under which the algorithm converges
% %in finitely many steps.



\begin{figure}[t]
  \small
    \centering
%initialization rule
\begin{prooftree}
\justifies
\volume \gets 0
\quad \Psi \gets \rect_\varphi
\using \hdecomp
\end{prooftree}

\vspace{.2in}
\begin{prooftree}
    m \models \Psi
    %\quad\quad
    %\emph{hvol} =  \int_{H^m} \prod_{x \in \vars_\varphi} \pdf(x)~ d\vars_\varphi
    \justifies
    \volume \gets \volume + \vol(H^m,\pdfs)
    \quad \Psi \gets \Psi \land \block(H^m)
    \using \sample
\end{prooftree}

\vspace{.1in}

$$\text{where } \block(H^m) \equiv \bigvee_{x \in \vars_\varphi}
  u_x < H^m_l(x) \lor l_x > H^m_u(x)$$
  \caption{\volalg:  weighted volume computation algorithm}
\label{alg:vol}
\end{figure}


\subsection{Weighted Volume Computation Algorithm}
\label{sec:swv}
To compute the integral over the region
$\varphi$, we exploit the observation
that if $\varphi$ is a \emph{hyperrectangular region},
i.e., an $n$-dimensional rectangle in $\mathds{R}^n$, 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 $\mathds{R}^2$:
$$\varphi \equiv 0 \leq x_1 \leq 100 \land 4 \leq x_2 \leq 10$$
%
\begin{flalign*}
    \text{The following holds:} ~~&\textstyle \int_\varphi \pdf_1(x_1)\pdf_2(x_2)~ dx_1dx_2 &\\
    =&\textstyle(\int_0^{100} \pdf_1(x_1)~ dx_1) (\int_4^{10} \pdf_2(x_2)~ dx_2) &\\
    =&\textstyle(\cdf_1(10) - \cdf_1(4))(\cdf_2(100) - \cdf_2(0))&
\end{flalign*}
%
where $\cdf_i$ is the \abr{CDF} of $\pdf_i(x_i)$,
which is $F_i(x) = \int^{x}_{-\infty} \pdf_i(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 \emph{hyperrectangular decomposition} phase
  represents the formula $\varphi$ 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 $\varphi$ symbolically as a formula $\rect_\varphi$, where each model
  of $\rect_\varphi$ corresponds to a hyperrectangle that
  lies inside the region $\varphi$.
%
  Second, after characterizing the set $\rect_\varphi$ of all
  hyperrectangles in $\varphi$, we can iteratively
  \emph{sample hyperrectangles} in $\varphi$,
  which can be done using an off-the-shelf \abr{SMT}
  solver to find models of $\rect_\varphi$.
  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 $\varphi$---that is, a lower bound on the exact
  weighted volume of $\varphi$.




\paragraph{Hyperrectangular decomposition}
We begin by defining hyperrectangles as  special formulas.
\begin{definition}[Hyperrectangles and their weighted volume]
  A formula $H \in \smt$ is a hyperrectangle
  if it can be written in the form
  $\bigwedge_{x \in \vars_H} c_x \leq x \leq c_x'$
  where $c_x,c'_x \in \mathds{R}$ are the lower
  and upper bounds of dimension $x$.
  We use $H_l(x)$ and $H_u(x)$
  to denote the lower and upper bounds of $x$ in $H$.

  The weighted volume of $H$, given a set $\pdfs$, is as follows:
  $$\textstyle \vol(H,\pdfs) = \prod_{x_i \in \vars_H} \int_{H_l(x_i)}^{H_u(x_i)} p_i(x_i) ~dx_i$$
\end{definition}



% \begin{figure}[t]
%   \centering
% \includegraphics[scale=1.2]{lowprob.pdf}
%
%
% \caption{Illustration of importance of directed sampling}
% \label{fig:decomp}
% \end{figure}


Ideally, we would take a formula $\varphi$
and rewrite it as a disjunction of hyperrectangles $\bigvee 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~\ref{fig:decomp}(a).
Here, there is no finite number of rectangles whose union is
the full region in $\mathds{R}^2$ enclosed by the triangle.

While the number of hyperrectangles enclosed in $\varphi$
is infinite, we can characterize them symbolically
using universal quantifiers, as shown by Li et al.~\cite{li2014management}.
Specifically, we define the hyperrectangular
decomposition of $\varphi$ as follows:

\begin{definition}[Hyperrectangular decomposition]
  \label{def:decomp}
  Given a formula $\varphi$, its \emph{hyperrectangular
  decomposition} $\rect_\varphi$ is:% defined as follows:
  $$\rect_\varphi \equiv \biggl(\bigwedge_{x \in \vars_\varphi} l_x \leq u_x\biggr) \land
    \forall \vars_\varphi \ldotp
      \biggl(\bigwedge_{x \in \vars_\varphi}
        l_x \leq x \leq u_x\biggr) \Rightarrow \varphi$$
  where $l_x,u_x$ are fresh free variables introduced for each $x \in \vars_\varphi$,
  and $\forall \vars_\varphi$ is short for $\forall x_1,\ldots,x_n$,
  for $x_i \in \vars_\varphi$.
%

Given a model $m \models \rect_\varphi$, we say
that $H^m$ is the \emph{hyperrectangle induced by} $m$,
as defined below:
$$H^m \equiv \bigwedge_{x\in \vars_\varphi} m(l_x) \leq x \leq m(u_x)$$
\end{definition}

Intuitively, $\rect_\varphi$ characterizes every possible
hyperrectangle that is subsumed by $\varphi$.
%Specifically, a model $m \models \rect_\varphi$ is
%over the set of free variables $\{l_{x_1}, u_{x_1},\ldots, l_{x_n},u_{x_n}\}$,
%where $\vars_\varphi = \{x_1, \ldots,x_n\}$.
The idea is that the hyperrectangle $H^m$ induced by each model
$m$ of $\rect_\varphi$ is %one of the hyperrectangles
%in $\varphi$.
%
%From the definition of $\rect_\varphi$, it is
%easy to see that $H^m$ is
subsumed by $\varphi$,
that is, $H^m \Rightarrow \varphi$.
%
The following example illustrates this process.
\begin{figure}[t]
  \centering
  \small
\includegraphics[scale=1]{triangle.pdf}

%
\hspace{0in}(a)\hspace{1.7in}(b)
%
\caption{$\mathds{R}^2$ view of (a) hyperrectangular decomposition
and (b) hyperrectangle sampling, where density is concentrated
in the top-left corner}
\label{fig:decomp}
\end{figure}





\begin{example}
  Consider the formula
  $\varphi \equiv x \geq y \land y \geq 0$,
  illustrated in Figure~\ref{fig:rectangle}
  as a gray, unbounded polyhedron.
  The formula $\rect_\varphi$, after eliminating the universal quantifier,
  is:
  $$l_x \leq u_x \land l_y \leq u_y \land l_y \geq 0 \land l_x \geq u_y$$
  Figure~\ref{fig:rectangle} shows two models $m_1,m_2 \models \rect_\varphi$
  and their graphical representation  as rectangles $H^{m_1},H^{m_2}$ in $\mathds{R}^2$. Observe that both rectangles are subsumed by $\varphi$.
\end{example}

% The following theorem states the soundness and completeness of hyperrectangular
% decomposition: models of $\rect_\varphi$ characterize all hyperrectangles in $\varphi$
% and no others.
%
% \begin{theorem}[Correctness of $\rect$]\label{thm:decomp}
% Let $\varphi \in \smt$.
% %\begin{itemize}
% %\item
%   \textbf{Soundness}:
% Let $m \models \rect_\varphi$. Then, $H^m \Rightarrow \varphi$
% is valid.
% %\item
% \textbf{Completeness}:
% Let $H$ be a hyperrectangle such that  $H \Rightarrow \varphi$.
% Then, the following is satisfiable:
% $\rect_\varphi \land \bigwedge_{x \in \vars_\varphi} l_x = H_l(x) \land u_x = H_u(x) $
% %\end{itemize}
% \end{theorem}

\paragraph{Hyperrectangle sampling}
Our symbolic weighted volume computation algorithm, $\volalg$,
is shown in Figure~\ref{alg:vol} as two
transition rules.
Given a pair $(\varphi, \pdfs)$,
the algorithm maintains a state
consisting of two variables: \rone $\volume$,
the current lower bound of the weighted
volume, and \rtwo $\Psi$, a constraint that encodes the \emph{remaining}
rectangles in the hyperrectangular decomposition of $\varphi$.

The algorithm is presented as guarded rules.
Initially, using the rule $\hdecomp$, $\volume$
is set to $0$ and $\Psi$ is set to $\rect_\varphi$.
%
The algorithm then proceeds by iteratively
applying the rule \sample.
Informally, the rule $\sample$ is used to
find arbitrary hyperrectangles in $\varphi$
and compute their weighted volume.
Specifically,
\sample finds a model $m$ of $\Psi$,
computes the weighted  volume of the hyperrectangle $H^m$ induced
by $m$, and adds the result to $\volume$.
%Recall that, by definition, $H^m \Rightarrow \varphi$.

To maintain soundness,
$\sample$ 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 constraint to $\Psi$---denoted $\block(H^m)$
and defined in Figure~\ref{alg:vol}---%
that ensures that for all models $m' \models\Psi$,
$H^{m'}$ does not  overlap with $H^m$, i.e., $H^{m'} \land H^{m}$
is unsatisfiable.
Informally, the $\block(H^m)$ constraint specifies
that any newly sampled hyperrectangle should be to
the \emph{left} or \emph{right} of $H^m$ for at least one of the dimensions.

% The following theorem  states the correctness of $\block$:
%  it removes all hyperrectangles that overlap with $H^m$
%  (soundness),
% and it does not overconstrain $\Psi$ by removing hyperrectangles
% that do not overlap with $H^m$ (completeness).
%
% \begin{theorem}[Correctness of $\block$] \label{thm:block}
% Given $\varphi$, let $\Psi \Rightarrow \rect_\varphi$, and
% let $m_1,m_2 \models \Psi$.
% %\begin{itemize}
% %
% %\item
% \textbf{Soundness}:
%     If $H^{m_1} \land H^{m_2}$ is satisfiable,
%     then $m_2 \not\models \Psi \land \block(H^{m_1})$.
% %
% %\item
%     \textbf{Completeness}: If $H^{m_1} \land H^{m_2}$
% is unsatisfiable,
% then $m_2 \models \Psi \land \block(H^{m_1})$.
% %\end{itemize}
% \end{theorem}

\begin{figure}[t]
    \centering
\includegraphics[scale=1]{decomp.pdf}

\caption{Illustration of models of $\rect_\varphi$}
\label{fig:rectangle}
\end{figure}

\paragraph{Lower and upper bounds}
The following theorem states the soundness of $\volalg$:
it maintains a lower bound on the exact
weighted volume.

\begin{theorem}[Soundness of $\volalg$]\label{thm:volsound}
The following is an invariant of
$\volalg(\varphi,\pdfs)$: $\volume \leq \vol(\varphi,\pdfs)$.
\end{theorem}

It follows from the above theorem that
we can use $\volalg$ to compute an upper bound
on the exact volume.
Specifically, because we are integrating over
\abrs{PDF}, we know that
$\vol(\varphi, \pdfs) + \vol(\neg \varphi,\pdfs) = 1$.
Therefore, by using \volalg to compute the weighted volume
of $\neg\varphi$, we get an upper bound
on the exact volume of $\varphi$.%, as stated in the following
%corollary.

\begin{corollary}[Computing upper bounds]
  \label{col:upper}
The following is an invariant of
$\volalg(\neg \varphi,\pdfs)$: $1 - \volume \geq \vol(\varphi,\pdfs)$
\end{corollary}



\subsection{Density-directed Sampling}
\label{sec:density}
While the \volalg algorithm is sound, it provides
no progress guarantees. Consider, for example,
a run that only samples unit hyperrectangles, i.e., points in $\mathds{R}^n$;
the volume of a point is 0, therefore, we will never compute
any volume.
Alternatively, the algorithm might diverge
by sampling hyperrectangles
in $\varphi$ that appear in very low probability density regions.
These two scenarios are illustrated in Figure~\ref{fig:decomp}(b)
on a triangular polytope in $\mathds{R}^2$.

Ideally, the rule \sample  would always find a model $m$ yielding
the hyperrectangle $H^m$ with the \emph{largest  weighted volume}.
Finding such a model amounts to solving the optimization
problem:
%
\begin{align*}
    \textstyle \underset{m \models \Psi}{\text{arg max}} \prod_{x_i \in \vars_\varphi}
    \int_{H^m_l(x_i)}^{H^m_u(x_i)}
    \pdf_i(x_i)~ dx_i
\end{align*}
%
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.

\begin{figure}[t]
  \centering
\begin{overpic}[scale=1.05]{gauss.pdf}
    \put(30,30) {\footnotesize (a)}
    \put(62,30) {\footnotesize (b)}
    \put(94,30) {\footnotesize (c)}
\end{overpic}
%\includegraphics[scale=1.05]{gauss.pdf}
%
%(a)\hspace{1in}(b)\hspace{.8in}(c)

\caption{Three \abrs{ADF} (gray) of a Gaussian \abr{PDF} (red) with mean 0 and
a standard deviation 3.
\abrs{ADF} are: (a) \emph{fine-grained};
(b) \emph{coarse}; (c) \emph{uniform}.}
\label{fig:approx}
\end{figure}

However, we make the key observation that
if $\pdf(x)$ is a \emph{step function}---i.e., piecewise constant---%
then we can symbolically encode the integral $\int \pdf(x)~dx$
in linear arithmetic.
As such, we propose to \rone \emph{approximate}
each density function $\pdf(x)$ with a step function $\appr(x)$,
\rtwo \emph{encode} the integrals $\int \appr(x)~dx$
as linear arithmetic formulas,
and \rthree \emph{direct} sampling towards hyperrectangles
that maximize these integrals, thus finding hyperrectangles
of large volume.%, as per the approximate density functions.

\paragraph{Approximate density functions and linear encodings}
%
We begin by defining \emph{approximate density functions} (\abrs{ADF}).
%as step functions.

\begin{definition}[Approximate density functions]
    An approximate density function  $\appr(x)$
    is of the following form:
 \[
    \appr(x) = \begin{cases}
        c_i, &  x \in [a_i,b_i) ~~~\text{for}~~~ 1 \leq i \leq n\\
        0, & \text{otherwise}
      \end{cases}
  \]
where %$n\in \mathds{N}$,
$c_i,a_i,b_i \in \mathds{R}$, $c_i > 0$,
and all $[a_i,b_i)$ are disjoint.
\end{definition}



We now show how to encode a formula $\sappr(x)$
over  the free variables $\delta_x, l_x, u_x$,
where for any model $m \models \sappr(x)$,
the value
 $m(\delta_x)$ is the area under $\appr(x)$
between $m(l_x)$ and $m(u_x)$, i.e.:
$m(\delta_x) = \int_{m(l_x)}^{m(u_x)} \appr(x)~dx$.
%
Intuitively, the value of this integral
is the sum of the areas of each \emph{bar} in $\appr(x)$,
restricted to%the interval
$[m(l_x), m(u_x)]$.

%\begin{definition}[Encoding area under an \abr{ADF}]
%Given an $\abr{ADF}$ $\appr(x)$, we define $\sappr(x)$ as follows:
%\[
%\sappr(x) \equiv \delta_x =  \sum_{i = 1}^n c_i* [a_i,b_i) \sqcap [l_x,u_x]
%\]
%where $\sqcap$ is the \emph{size of the intersection} of two intervals:
%\[
%[a_i,b_i) \sqcap [l_x,u_x] = \text{max}(\text{min}(b_i,u_x) - \text{max}(a_i,l_x),0)
%\]
%\end{definition}

\begin{definition}[Encoding area under an \abr{ADF}]
Given an $\abr{ADF}$ $\appr(x)$, we define $\sappr(x)$ as follows:
\[
    \textstyle \sappr(x) \equiv \delta_x = \sum_{i = 1}^n c_i \cdot \bigl|[a_i,b_i) \cap [l_x,u_x]\bigr|
\]
\end{definition}

The  constraint $\sappr(x)$ is directly expressible
in linear arithmetic,
since
\[
\bigl|[a_i,b_i) \cap [l_x,u_x]\bigr| = \text{max}(\text{min}(b_i,u_x) - \text{max}(a_i,l_x),0)
\]
The finite sum in $\sappr(x)$ computes the size
of the intersection of $[l_x,u_x]$ with each
interval $[a_i,b_i)$ in  $\appr(x)$, and multiplies
the intersection with $c_i$, the value of the $\appr$ in that interval.
% The following theorem states the correctness of the above
% encoding:
%
% \begin{theorem}[Correctness of $\sappr$]\label{thm:adf}
% Fix an \abr{ADF} $\appr(x)$.
% %\begin{itemize}
%   %\item
%   \textbf{Soundness}:
% For any model $m \models \sappr$,
% the following  is true:
% $m(\delta_x) = \int_{m(l_x)}^{m(u_x)} \appr(x)~dx$.
% %\item
% \textbf{Completeness}:
% For any constants $a,b,c\in \mathds{R}$ such that
% $c = \int_a^b \appr(x)~dx$, the following formula is satisfiable:
% $\delta_x = c \land l_x = a \land u_x = b \land \sappr(x)$.
% %\end{itemize}
% \end{theorem}





\begin{figure}[t!]
  \footnotesize
    \centering
%initialization rule
\begin{prooftree}
\justifies
\volume \gets 0
\quad \Psi \gets \rect_\varphi
\quad \emph{lb} \gets 1
\using \hdecomp
\end{prooftree}

\vspace{.2in}
\begin{prooftree}
    m \models \Psi \land \bigwedge_{x_i\in \vars_\varphi}\exists\delta_{x_i} \ldotp \sappr_i(x_i) \land \delta_{x_i} \geq \emph{lb}
    %\quad
    %\emph{hvol} =  \int_{H^m} \prod_{x \in \vars_\varphi} \pdf(x)~ %d\vars_\varphi
    \justifies
    \volume \gets \volume +\vol(H^m,\pdfs)
    \quad\quad \Psi \gets \Psi \land \block(H^m)
    \using \sample
\end{prooftree}

\vspace{.1in}
\begin{prooftree}
    \justifies
    \emph{lb} \gets \lambda*\emph{lb}
    \using \decay
\end{prooftree}


\caption{\volopt: \abr{ADF}-directed volume computation}
\label{fig:algadf}
\end{figure}

\paragraph{\abr{ADF}-directed Volume Computation}
We now present the algorithm \volopt (Figure~\ref{fig:algadf}), an extension of our volume computation
algorithm \volalg that uses \abrs{ADF} to steer the sampling process.
The \abrs{ADF} are only used for guiding the rule
\sample towards dense hyperrectangles, and thus do not affect
soundness of the volume computation.
For example, Figure~\ref{fig:approx} shows three approximations
of a Gaussian; all three are valid approximations.
In Section~\ref{sec:impl}, we discuss the impact
of different \abrs{ADF} on performance.

Formally, we create a set of \abrs{ADF} $\apprs = \{\appr_1,\ldots,\appr_n\}$,
where, for each variable $x_i \in \vars_\varphi$,
we  associate the \abr{ADF} $\appr_i(x_i)$.
The rule \sample now encodes $\sappr_i(x_i)$ and attempts
to find a hyperrectangle such that for each dimension $x$,
$\delta_x$ is greater than some lower bound $\emph{lb}$,
which is initialized to $1$.
Of course, we need to keep reducing the value $\emph{lb}$ as we run
out of hyperrectangles of a given volume.
Therefore, the rule \decay is used to shrink $\emph{lb}$ using a
fixed \emph{decay rate} $\lambda \in (0,1)$
and can be applied when \sample fails to find
a sufficiently large hyperrectangle.

Note that, ideally, we would look for a model $m$
such that $\prod_{x\in \vars_\varphi} \delta_x$ is maximized,
thus, finding the hyperrectangle with the largest weighted volume
with respect to the \abrs{ADF}.
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 $\delta_x$
is greater than the lower bound.





\subsection{Convergence of \volopt}
We now discuss the convergence properties of $\volopt$.
Suppose we are given a  formula  $\varphi$, a set $\pdfs$,
and a set $\apprs$.
Let $R \subset \mathds{R}^n$ be the region where all the \abrs{ADF} in $\apprs$
are non-zero.
We will show that $\volopt$ converges, in the limit, to the
exact weighted volume restricted to $R$;
that is, $\volopt$ converges to
%
$\textstyle \int_{\varphi \cap R} \prod_{x_i\in \vars_\varphi} \pdf_i(x_i)~ d\vars_\varphi$.

Note the fascinating part here is that we do not impose
any restrictions on the \abrs{ADF}:
they do not have to
have any correspondence with the \abrs{PDF}
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~\ref{sec:impl}.

The following theorem states convergence of
$\volopt$; it assumes that \sample is applied
iteratively and \decay is only
applied when \sample cannot find a model.

\begin{theorem}[Convergence to $R$]\label{thm:limit}
    Assume $\volopt$ is run on $(\varphi,\pdfs)$
    and a set of \abrs{ADF}  $\apprs$
that are non-zero for $R \subset \mathds{R}^n$.
Let $\volume_i$ be the value of $\volume$ after
$i$ applications of \sample. Then,
$\textstyle \lim_{i \rightarrow \infty} \volume_i =  \int_{\varphi \cap R} \prod_{x_i\in \vars_\varphi} \pdf_i(x_i)~ d\vars_\varphi$.
\end{theorem}

% A full proof of this theorem may be found in the \iftoggle{long}{Appendix}{\emph{supplementary materials}};
% here we provide a sketch of the intuition:
% suppose, for the sake of obtaining a contradiction,
% that we never sample some hyperrectangle $H$ inside of $R$.
% $H$ has a non-zero weighted volume with respect to the approximation,
% so we \emph{would} sample $H$ before finitely many applications of \decay.
% Thus, we sample all of $R$ in the limit, and
% $\volume_i \rightarrow \int_{\varphi \cap R} \prod_{x_i \in \vars_\varphi} \pdf_i(x_i)~ d\vars_\varphi$.
%
Note that the above theorem
directly gives us a way to approach the exact volume.
Specifically, by performing runs of $\volopt$
on subsets in an infinite partition of $\mathds{R}^n$
induced by the \abrs{ADF},
we can ensure that the sum over the
$\volopt$ processes approaches the exact volume in the limit.
For all $i$, let $\apprs_i$ be a set of \abrs{ADF}
corresponding to an $\volopt$ process $P_i$,
where $R_i \subset \mathds{R}^n$ is the non-zero
region of $\apprs_i$.
We require an infinite set of $P_i$ to partition $\mathds{R}^n$:
\rone for all $i \neq j$, $R_i \cap R_j = \emptyset$, and
\rtwo $\bigcup_{i=1}^\infty R_i = \mathds{R}^n$.
The following theorem formalizes the argument:

\begin{theorem}[Convergence]\label{thm:limitfull}
  Let $P_1,P_2,\ldots$ be $\volopt$ processes
  that partition $\mathds{R}^n$.
  Assume a fair serialization
  where each $P_i$ performs $\sample$
  infinitely often, and
  let $\volume_n$ be the total computed volume
  across all $P_i$ after $n$ successful calls to $\sample$.
  Then, $\lim_{n \rightarrow \infty} \volume_n = \vol(\varphi,\pdfs)$.
\end{theorem}




%!TEX root = paper.tex

paper.bbl
\begin{thebibliography}{72}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[dat()]{dataset}
Uci machine learning repository: Census income.
\newblock \url{https://archive.ics.uci.edu/ml/datasets/Adult/}.

\bibitem[mat({\natexlab{a}})]{mathematica}
Wolfram mathematica.
\newblock \url{https://www.wolfram.com/mathematica/}, {\natexlab{a}}.

\bibitem[mat({\natexlab{b}})]{matlab}
Matlab.
\newblock \url{http://www.mathworks.com/products/matlab/}, {\natexlab{b}}.

\bibitem[red()]{redlog}
Redlog.
\newblock \url{http://www.redlog.eu/}.

\bibitem[wek()]{weka}
Weka.
\newblock \url{http://www.cs.waikato.ac.nz/ml/weka/}.

\bibitem[EEO(2014)]{EEOC}
Code of federal regulations.
\newblock
  \url{https://www.gpo.gov/fdsys/pkg/CFR-2014-title29-vol4/xml/CFR-2014-title29-vol4-part1607.xml},
  July 2014.
\newblock (Accessed on 06/18/2016).

\bibitem[Ajunwa et~al.(2016)Ajunwa, Friedler, Scheidegger, and
  Venkatasubramanian]{ajunwa2016hiring}
I.~Ajunwa, S.~Friedler, C.~E. Scheidegger, and S.~Venkatasubramanian.
\newblock Hiring by algorithm: predicting and preventing disparate impact.
\newblock \emph{Available at SSRN 2746078}, 2016.

\bibitem[Angwin et~al.(2016)Angwin, Larson, Mattu, and
  Kirchner]{articlePropublica}
J.~Angwin, J.~Larson, S.~Mattu, and L.~Kirchner.
\newblock Machine bias: There's software used across the country to predict
  future criminals. and it's biased against blacks.
\newblock
  \url{https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-sentencing},
  May 2016.
\newblock (Accessed on 06/18/2016).

\bibitem[Asarin et~al.(2000)Asarin, Bournez, Dang, and
  Maler]{asarin2000approximate}
E.~Asarin, O.~Bournez, T.~Dang, and O.~Maler.
\newblock Approximate reachability analysis of piecewise-linear dynamical
  systems.
\newblock In \emph{International Workshop on Hybrid Systems: Computation and
  Control}, pages 20--31. Springer, 2000.

\bibitem[Attali and Burstein(2006)]{attali06}
Y.~Attali and J.~Burstein.
\newblock Automated essary scoring with e-rater v.2.
\newblock \emph{Journal of Technology, Learning, and Assessment}, 4\penalty0
  (3), 2006.

\bibitem[Barnett and Leino(2005)]{barnett2005weakest}
M.~Barnett and K.~R.~M. Leino.
\newblock Weakest-precondition of unstructured programs.
\newblock In \emph{ACM SIGSOFT Software Engineering Notes}, volume~31, pages
  82--87. ACM, 2005.

\bibitem[Barocas and Selbst(2014)]{barocas2014big}
S.~Barocas and A.~D. Selbst.
\newblock Big data's disparate impact.
\newblock \emph{Available at SSRN 2477899}, 2014.

\bibitem[Barthe et~al.(2004)Barthe, D'Argenio, and Rezk]{barthe04}
G.~Barthe, P.~R. D'Argenio, and T.~Rezk.
\newblock Secure information flow by self-composition.
\newblock In \emph{CSFW}, 2004.

\bibitem[Barthe et~al.(2014)Barthe, Gaboardi, Arias, Hsu, Kunz, and
  Strub]{barthe14}
G.~Barthe, M.~Gaboardi, E.~J.~G. Arias, J.~Hsu, C.~Kunz, and P.~Strub.
\newblock Proving differential privacy in hoare logic.
\newblock In \emph{{IEEE} 27th Computer Security Foundations Symposium, {CSF}
  2014, Vienna, Austria, 19-22 July, 2014}, pages 411--424, 2014.
\newblock \doi{10.1109/CSF.2014.36}.
\newblock URL \url{http://dx.doi.org/10.1109/CSF.2014.36}.

\bibitem[Bastani et~al.(2016)Bastani, Ioannou, Lampropoulos, Vytiniotis, Nori,
  and Criminisi]{bastani16}
O.~Bastani, Y.~Ioannou, L.~Lampropoulos, D.~Vytiniotis, A.~Nori, and
  A.~Criminisi.
\newblock Measuring neural net robustness with constraints.
\newblock \emph{CoRR}, abs/1605.07262, 2016.
\newblock URL \url{http://arxiv.org/abs/1605.07262}.

\bibitem[Belle et~al.(2015{\natexlab{a}})Belle, Passerini, and den
  Broeck]{belle15}
V.~Belle, A.~Passerini, and G.~V. den Broeck.
\newblock Probabilistic inference in hybrid domains by weighted model
  integration.
\newblock In \emph{Proceedings of the Twenty-Fourth International Joint
  Conference on Artificial Intelligence, {IJCAI} 2015, Buenos Aires, Argentina,
  July 25-31, 2015}, pages 2770--2776, 2015{\natexlab{a}}.
\newblock URL \url{http://ijcai.org/Abstract/15/392}.

\bibitem[Belle et~al.(2015{\natexlab{b}})Belle, Van~den Broeck, and
  Passerini]{belle2015hashing}
V.~Belle, G.~Van~den Broeck, and A.~Passerini.
\newblock Hashing-based approximate probabilistic inference in hybrid domains.
\newblock In \emph{Proceedings of the 31st Conference on Uncertainty in
  Artificial Intelligence (UAI)}, 2015{\natexlab{b}}.

\bibitem[Belle et~al.(2016)Belle, den Broeck, and Passerini]{belle16}
V.~Belle, G.~V. den Broeck, and A.~Passerini.
\newblock Component caching in hybrid domains with piecewise polynomial
  densities.
\newblock In \emph{Proceedings of the Thirtieth {AAAI} Conference on Artificial
  Intelligence, February 12-17, 2016, Phoenix, Arizona, {USA.}}, pages
  3369--3375, 2016.
\newblock URL
  \url{http://www.aaai.org/ocs/index.php/AAAI/AAAI16/paper/view/12542}.

\bibitem[Berg(2014)]{articleGuardianCrime}
N.~Berg.
\newblock Predicting crime, lapd-style.
\newblock \url{https://www.theguardian.com/cities/2014/jun/25/predicting-
  crime-lapd-los-angeles-police-data-analysis-algorithm-minority-report}, June
  2014.
\newblock (Accessed on 06/18/2016).

\bibitem[Bishop(2006)]{bishop2006pattern}
C.~M. Bishop.
\newblock Pattern recognition.
\newblock \emph{Machine Learning}, 128, 2006.

\bibitem[Bournez et~al.(1999)Bournez, Maler, and Pnueli]{Bournez99}
O.~Bournez, O.~Maler, and A.~Pnueli.
\newblock \emph{Orthogonal Polyhedra: Representation and Computation}, pages
  46--60.
\newblock Springer Berlin Heidelberg, Berlin, Heidelberg, 1999.
\newblock ISBN 978-3-540-48983-2.
\newblock \doi{10.1007/3-540-48983-5_8}.
\newblock URL \url{http://dx.doi.org/10.1007/3-540-48983-5_8}.

\bibitem[Calders and Verwer(2010)]{calders10}
T.~Calders and S.~Verwer.
\newblock Three naive bayes approaches for discrimination-free classification.
\newblock \emph{Data Mining and Knowledge Discovery}, 21\penalty0 (2):\penalty0
  277--292, 2010.

\bibitem[Carbin et~al.(2013)Carbin, Misailovic, and Rinard]{carbin13}
M.~Carbin, S.~Misailovic, and M.~C. Rinard.
\newblock Verifying quantitative reliability for programs that execute on
  unreliable hardware.
\newblock In \emph{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.
\newblock \doi{10.1145/2509136.2509546}.
\newblock URL \url{http://doi.acm.org/10.1145/2509136.2509546}.

\bibitem[Chaudhuri et~al.(2011)Chaudhuri, Gulwani, Lublinerman, and
  Navidpour]{Chaudhuri11}
S.~Chaudhuri, S.~Gulwani, R.~Lublinerman, and S.~Navidpour.
\newblock Proving programs robust.
\newblock In \emph{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.
\newblock ISBN 978-1-4503-0443-6.
\newblock \doi{10.1145/2025113.2025131}.
\newblock URL \url{http://doi.acm.org/10.1145/2025113.2025131}.

\bibitem[Chaudhuri et~al.(2014)Chaudhuri, Clochard, and
  Solar-Lezama]{chaudhuri14}
S.~Chaudhuri, M.~Clochard, and A.~Solar-Lezama.
\newblock Bridging boolean and quantitative synthesis using smoothed proof
  search.
\newblock In \emph{POPL}, volume~49, pages 207--220. ACM, 2014.

\bibitem[Chib and Greenberg(1995)]{chib1995understanding}
S.~Chib and E.~Greenberg.
\newblock Understanding the metropolis-hastings algorithm.
\newblock \emph{The american statistician}, 49\penalty0 (4):\penalty0 327--335,
  1995.

\bibitem[Chistikov et~al.(2014)Chistikov, Dimitrova, and
  Majumdar]{chistikov14arxiv}
D.~Chistikov, R.~Dimitrova, and R.~Majumdar.
\newblock Approximate counting in {SMT} and value estimation for probabilistic
  programs.
\newblock \emph{CoRR}, abs/1411.0659, 2014.
\newblock URL \url{http://arxiv.org/abs/1411.0659}.

\bibitem[Chistikov et~al.(2015)Chistikov, Dimitrova, and Majumdar]{chistikov15}
D.~Chistikov, R.~Dimitrova, and R.~Majumdar.
\newblock Approximate counting in {SMT} and value estimation for probabilistic
  programs.
\newblock In \emph{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.
\newblock \doi{10.1007/978-3-662-46681-0_26}.
\newblock URL \url{http://dx.doi.org/10.1007/978-3-662-46681-0_26}.

\bibitem[Claret et~al.(2013)Claret, Rajamani, Nori, Gordon, and
  Borgstr{\"o}m]{claret2013bayesian}
G.~Claret, S.~K. Rajamani, A.~V. Nori, A.~D. Gordon, and J.~Borgstr{\"o}m.
\newblock Bayesian inference using data flow analysis.
\newblock In \emph{Proceedings of the 2013 9th Joint Meeting on Foundations of
  Software Engineering}, pages 92--102. ACM, 2013.

\bibitem[Clarke et~al.(2004)Clarke, Kroening, and Lerda]{clarke2004tool}
E.~Clarke, D.~Kroening, and F.~Lerda.
\newblock A tool for checking ansi-c programs.
\newblock In \emph{International Conference on Tools and Algorithms for the
  Construction and Analysis of Systems}, pages 168--176. Springer, 2004.

\bibitem[Cytron et~al.(1991)Cytron, Ferrante, Rosen, Wegman, and
  Zadeck]{cytron1991efficiently}
R.~Cytron, J.~Ferrante, B.~K. Rosen, M.~N. Wegman, and F.~K. Zadeck.
\newblock Efficiently computing static single assignment form and the control
  dependence graph.
\newblock \emph{ACM Transactions on Programming Languages and Systems
  (TOPLAS)}, 13\penalty0 (4):\penalty0 451--490, 1991.

\bibitem[Datta et~al.(2015)Datta, Tschantz, and Datta]{datta2015automated}
A.~Datta, M.~C. Tschantz, and A.~Datta.
\newblock Automated experiments on ad privacy settings.
\newblock \emph{Proceedings on Privacy Enhancing Technologies}, 2015\penalty0
  (1):\penalty0 92--112, 2015.

\bibitem[Datta et~al.(2016)Datta, Sen, and Zick]{datta2016algorithmic}
A.~Datta, S.~Sen, and Y.~Zick.
\newblock Algorithmic transparency via quantitative input influence.
\newblock In \emph{Proceedings of 37th IEEE Symposium on Security and Privacy},
  2016.

\bibitem[De~Loera et~al.(2012)De~Loera, Dutra, Koeppe, Moreinis, Pinto, and
  Wu]{de2012software}
J.~De~Loera, B.~Dutra, M.~Koeppe, S.~Moreinis, G.~Pinto, and J.~Wu.
\newblock Software for exact integration of polynomials over polyhedra.
\newblock \emph{ACM Communications in Computer Algebra}, 45\penalty0
  (3/4):\penalty0 169--172, 2012.

\bibitem[De~Moura and Bj{\o}rner(2008)]{de2008z3}
L.~De~Moura and N.~Bj{\o}rner.
\newblock Z3: An efficient smt solver.
\newblock In \emph{International conference on Tools and Algorithms for the
  Construction and Analysis of Systems}, pages 337--340. Springer, 2008.

\bibitem[Dwork(2006)]{dwork2006differential}
C.~Dwork.
\newblock Differential privacy.
\newblock In \emph{Automata, languages and programming}, pages 1--12. Springer,
  2006.

\bibitem[Dwork et~al.(2012)Dwork, Hardt, Pitassi, Reingold, and Zemel]{dwork12}
C.~Dwork, M.~Hardt, T.~Pitassi, O.~Reingold, and R.~S. Zemel.
\newblock Fairness through awareness.
\newblock In \emph{Innovations in Theoretical Computer Science 2012, Cambridge,
  MA, USA, January 8-10, 2012}, pages 214--226, 2012.
\newblock \doi{10.1145/2090236.2090255}.

\bibitem[Dyer et~al.(1991)Dyer, Frieze, and Kannan]{dyer1991random}
M.~Dyer, A.~Frieze, and R.~Kannan.
\newblock A random polynomial-time algorithm for approximating the volume of
  convex bodies.
\newblock \emph{Journal of the ACM (JACM)}, 38\penalty0 (1):\penalty0 1--17,
  1991.

\bibitem[Dyer and Frieze(1988)]{dyer1988complexity}
M.~E. Dyer and A.~M. Frieze.
\newblock On the complexity of computing the volume of a polyhedron.
\newblock \emph{SIAM Journal on Computing}, 17\penalty0 (5):\penalty0 967--974,
  1988.

\bibitem[Eubanks(2015)]{articleSlateWelfare}
V.~Eubanks.
\newblock The dangers of letting algorithms enforce policy.
\newblock
  \url{http://www.slate.com/articles/technology/future_tense/2015/04/the_dangers_of_letting_algorithms_enforce_policy.html},
  April 2015.
\newblock (Accessed on 06/18/2016).

\bibitem[Feldman et~al.(2015)Feldman, Friedler, Moeller, Scheidegger, and
  Venkatasubramanian]{feldman15}
M.~Feldman, S.~A. Friedler, J.~Moeller, C.~Scheidegger, and
  S.~Venkatasubramanian.
\newblock Certifying and removing disparate impact.
\newblock In \emph{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.
\newblock \doi{10.1145/2783258.2783311}.
\newblock URL \url{http://doi.acm.org/10.1145/2783258.2783311}.

\bibitem[Filieri et~al.(2013)Filieri, P{\u{a}}s{\u{a}}reanu, and
  Visser]{filieri2013reliability}
A.~Filieri, C.~S. P{\u{a}}s{\u{a}}reanu, and W.~Visser.
\newblock Reliability analysis in symbolic pathfinder.
\newblock In \emph{Proceedings of the 2013 International Conference on Software
  Engineering}, pages 622--631. IEEE Press, 2013.

\bibitem[Friedler et~al.(2016)Friedler, Scheidegger, and
  Venkatasubramanian]{friedler16}
S.~A. Friedler, C.~Scheidegger, and S.~Venkatasubramanian.
\newblock On the (im)possibility of fairness.
\newblock \emph{CoRR}, abs/1609.07236, 2016.
\newblock URL \url{http://arxiv.org/abs/1609.07236}.

\bibitem[Gehr et~al.(2016)Gehr, Misailovic, and Vechev]{Gehr16}
T.~Gehr, S.~Misailovic, and M.~Vechev.
\newblock Psi: Exact symbolic inference for probabilistic programs.
\newblock In \emph{Computer aided verification}. Springer, 2016.

\bibitem[Geldenhuys et~al.(2012)Geldenhuys, Dwyer, and
  Visser]{geldenhuys2012probabilistic}
J.~Geldenhuys, M.~B. Dwyer, and W.~Visser.
\newblock Probabilistic symbolic execution.
\newblock In \emph{Proceedings of the 2012 International Symposium on Software
  Testing and Analysis}, pages 166--176. ACM, 2012.

\bibitem[Goodman et~al.(2008)Goodman, Mansinghka, Roy, Bonawitz, and
  Tenenbaum]{goodman08}
N.~D. Goodman, V.~K. Mansinghka, D.~M. Roy, K.~Bonawitz, and J.~B. Tenenbaum.
\newblock Church: a language for generative models.
\newblock In \emph{{UAI} 2008, Proceedings of the 24th Conference in
  Uncertainty in Artificial Intelligence, Helsinki, Finland, July 9-12, 2008},
  pages 220--229, 2008.

\bibitem[Gordon et~al.(2014)Gordon, Henzinger, Nori, and
  Rajamani]{gordon2014probabilistic}
A.~D. Gordon, T.~A. Henzinger, A.~V. Nori, and S.~K. Rajamani.
\newblock Probabilistic programming.
\newblock In \emph{Proceedings of the on Future of Software Engineering}, pages
  167--181. ACM, 2014.

\bibitem[Hardt et~al.(2016)Hardt, Price, and Srebro]{hardt16}
M.~Hardt, E.~Price, and N.~Srebro.
\newblock Equality of opportunity in supervised learning.
\newblock \emph{CoRR}, abs/1610.02413, 2016.
\newblock URL \url{http://arxiv.org/abs/1610.02413}.

\bibitem[Khachiyan(1993)]{khachiyan1993complexity}
L.~Khachiyan.
\newblock \emph{Complexity of polytope volume computation}.
\newblock Springer, 1993.

\bibitem[Kobie(2016)]{articleWired}
N.~Kobie.
\newblock Who do you blame when an algorithm gets you fired?
\newblock \url{http://www.wired.co.uk/article/make-algorithms-accountable},
  January 2016.
\newblock (Accessed on 06/18/2016).

\bibitem[Koller and Friedman(2009)]{koller2009probabilistic}
D.~Koller and N.~Friedman.
\newblock \emph{Probabilistic graphical models: principles and techniques}.
\newblock MIT press, 2009.

\bibitem[Kozen(1981)]{kozen1981semantics}
D.~Kozen.
\newblock Semantics of probabilistic programs.
\newblock \emph{Journal of Computer and System Sciences}, 22\penalty0
  (3):\penalty0 328--350, 1981.

\bibitem[Li et~al.(2014)Li, Tan, and Chechik]{li2014management}
Y.~Li, T.~H. Tan, and M.~Chechik.
\newblock Management of time requirements in component-based systems.
\newblock In \emph{FM 2014: Formal Methods}, pages 399--415. Springer, 2014.

\bibitem[Mardziel et~al.(2011)Mardziel, Magill, Hicks, and
  Srivatsa]{mardziel2011dynamic}
P.~Mardziel, S.~Magill, M.~Hicks, and M.~Srivatsa.
\newblock Dynamic enforcement of knowledge-based security policies.
\newblock In \emph{Computer Security Foundations Symposium (CSF), 2011 IEEE
  24th}, pages 114--128. IEEE, 2011.

\bibitem[Miller(2015)]{articleTimesHiring}
C.~C. Miller.
\newblock Can an algorithm hire better than a human?
\newblock
  \url{http://www.nytimes.com/2015/06/26/upshot/can-an-algorithm-hire-better-than-a-human.html},
  June 2015.
\newblock (Accessed on 06/18/2016).

\bibitem[Minka et~al.(2012)Minka, Winn, Guiver, and Knowles]{minka2012infer}
T.~Minka, J.~Winn, J.~Guiver, and D.~Knowles.
\newblock Infer.net 2.5.
\newblock \emph{Microsoft Research Cambridge}, 2012.

\bibitem[Monniaux(2000)]{monniaux2000abstract}
D.~Monniaux.
\newblock Abstract interpretation of probabilistic semantics.
\newblock In \emph{Static Analysis}, pages 322--339. Springer, 2000.

\bibitem[Monniaux(2001{\natexlab{a}})]{monniaux2001abstract}
D.~Monniaux.
\newblock An abstract monte-carlo method for the analysis of probabilistic
  programs.
\newblock In \emph{ACM SIGPLAN Notices}, volume~36, pages 93--101. ACM,
  2001{\natexlab{a}}.

\bibitem[Monniaux(2001{\natexlab{b}})]{monniaux2001backwards}
D.~Monniaux.
\newblock Backwards abstract interpretation of probabilistic programs.
\newblock In \emph{Programming Languages and Systems}, pages 367--382.
  Springer, 2001{\natexlab{b}}.

\bibitem[Nair and Hinton(2010)]{nair2010rectified}
V.~Nair and G.~E. Hinton.
\newblock Rectified linear units improve restricted boltzmann machines.
\newblock In \emph{Proceedings of the 27th International Conference on Machine
  Learning (ICML-10)}, pages 807--814, 2010.

\bibitem[Pedreshi et~al.(2008)Pedreshi, Ruggieri, and
  Turini]{pedreshi2008discrimination}
D.~Pedreshi, S.~Ruggieri, and F.~Turini.
\newblock Discrimination-aware data mining.
\newblock In \emph{Proceedings of the 14th ACM SIGKDD international conference
  on Knowledge discovery and data mining}, pages 560--568. ACM, 2008.

\bibitem[Perry(2013)]{perry2013predictive}
W.~L. Perry.
\newblock \emph{Predictive policing: The role of crime forecasting in law
  enforcement operations}.
\newblock Rand Corporation, 2013.

\bibitem[Rawls(2009)]{rawls2009theory}
J.~Rawls.
\newblock \emph{A theory of justice}.
\newblock Harvard university press, 2009.

\bibitem[Ruggieri(2014)]{ruggieri2014using}
S.~Ruggieri.
\newblock Using t-closeness anonymity to control for non-discrimination.
\newblock \emph{Transactions on Data Privacy}, 7\penalty0 (2):\penalty0
  99--129, 2014.

\bibitem[Sampson et~al.(2014)Sampson, Panchekha, Mytkowicz, McKinley, Grossman,
  and Ceze]{sampson2014expressing}
A.~Sampson, P.~Panchekha, T.~Mytkowicz, K.~S. McKinley, D.~Grossman, and
  L.~Ceze.
\newblock Expressing and verifying probabilistic assertions.
\newblock In \emph{ACM SIGPLAN Notices}, volume~49, pages 112--122. ACM, 2014.

\bibitem[Sankaranarayanan et~al.(2013)Sankaranarayanan, Chakarov, and
  Gulwani]{sankaranarayanan13}
S.~Sankaranarayanan, A.~Chakarov, and S.~Gulwani.
\newblock Static analysis for probabilistic programs: inferring whole program
  properties from finitely many paths.
\newblock In \emph{{ACM} {SIGPLAN} Conference on Programming Language Design
  and Implementation, {PLDI} '13, Seattle, WA, USA, June 16-19, 2013}, pages
  447--458, 2013.
\newblock \doi{10.1145/2462156.2462179}.
\newblock URL \url{http://doi.acm.org/10.1145/2462156.2462179}.

\bibitem[Sweeney(2013)]{sweeney2013discrimination}
L.~Sweeney.
\newblock Discrimination in online ad delivery.
\newblock \emph{Queue}, 11\penalty0 (3):\penalty0 10, 2013.

\bibitem[Tutt(2016)]{tutt2016fda}
A.~Tutt.
\newblock An fda for algorithms.
\newblock \emph{Available at SSRN 2747994}, 2016.

\bibitem[Valentino-Devries et~al.(2012)Valentino-Devries, Singer-Vine, and
  Soltani]{articWsjstaples}
J.~Valentino-Devries, J.~Singer-Vine, and A.~Soltani.
\newblock Websites vary prices, deals based on users' information.
\newblock
  \url{http://www.wsj.com/articles/SB10001424127887323777204578189391813881534},
  December 2012.
\newblock (Accessed on 06/18/2016).

\bibitem[Vempala(2005)]{vempala2005geometric}
S.~Vempala.
\newblock Geometric random walks: a survey.
\newblock \emph{Combinatorial and computational geometry}, 52\penalty0
  (573-612):\penalty0 2, 2005.

\bibitem[Zarsky(2014)]{zarsky2014understanding}
T.~Zarsky.
\newblock Understanding discrimination in the scored society.
\newblock \emph{Washington Law Review}, 89\penalty0 (4), 2014.

\bibitem[Zemel et~al.(2013)Zemel, Wu, Swersky, Pitassi, and Dwork]{zemel13}
R.~S. Zemel, Y.~Wu, K.~Swersky, T.~Pitassi, and C.~Dwork.
\newblock Learning fair representations.
\newblock In \emph{Proceedings of the 30th International Conference on Machine
  Learning, {ICML} 2013, Atlanta, GA, USA, 16-21 June 2013}, pages 325--333,
  2013.
\newblock URL \url{http://jmlr.org/proceedings/papers/v28/zemel13.html}.

\end{thebibliography}
