# Solving Satisfiability Modulo Counting for Symbolic and Statistical AI
Integration With Provable Guarantees

Jinzhao Li, Nan Jiang, Yexiang Xue
Department of Computer Science, Purdue University
<EMAIL_ADDRESS>
###### Abstract

Satisfiability Modulo Counting (SMC) encompasses problems that require both
symbolic decision-making and statistical reasoning. Its general formulation
captures many real-world problems at the intersection of symbolic and
statistical Artificial Intelligence. SMC searches for policy interventions to
control probabilistic outcomes. Solving SMC is challenging because of its
highly intractable nature ($\text{NP}^{\text{PP}}$-complete), incorporating
statistical inference and symbolic reasoning. Previous research on SMC solving
lacks provable guarantees and/or suffers from sub-optimal empirical
performance, especially when combinatorial constraints are present. We propose
$\mathtt{XOR}$-$\mathtt{SMC}$, a polynomial algorithm with access to NP-
oracles, to solve highly intractable SMC problems with constant approximation
guarantees. $\mathtt{XOR}$-$\mathtt{SMC}$ transforms the highly intractable
SMC into satisfiability problems, by replacing the model counting in SMC with
SAT formulae subject to randomized XOR constraints. Experiments on solving
important SMC problems in AI for social good demonstrate that
$\mathtt{XOR}$-$\mathtt{SMC}$ finds solutions close to the true optimum,
outperforming several baselines which struggle to find good approximations for
the intractable model counting in SMC.

##  1 Introduction

Symbolic and statistical approaches are two fundamental driving forces of
Artificial Intelligence (AI). Symbolic AI, exemplified by SATisfiability (SAT)
and constraint programming, finds solutions satisfying constraints but
requires rigid formulations and is difficult to include probabilities.
Statistical AI captures uncertainty but often lacks constraint satisfaction.
Integrating symbolic and statistical AI remains an open field and has gained
research attention recently [1, 2, 3].

Satisfiability Modulo Counting (SMC) is an umbrella problem at the
intersection of symbolic and statistical AI. It encompasses problems that
carry out symbolic decision-making (satisfiability) mixed with statistical
reasoning (model counting). SMC searches for policy interventions to control
probabilistic outcomes. Formally, SMC is an SAT problem involving predicates
on model counts. Model counting computes the number of models (i.e.,
solutions) to an SAT formula. Its weighted form subsumes probabilistic
inference on Machine Learning (ML) models.

As a motivating SMC application, stochastic connectivity optimization searches
for the optimal plan to reinforce the network structure so its connectivity
preserves under stochastic events – a central problem for a city planner who
works on securing her residents multiple paths to emergency shelters in case
of natural disasters. This problem is useful for disaster preparation [4],
bio-diversity protection [5], internet resilience [6], social influence
maximization [7], energy security [8], etc. It requires symbolic reasoning
(satisfiability) to decide which roads to reinforce and where to place
emergency shelters, and statistical inference (model counting) to reason about
the number of paths to shelters and the probabilities of natural disasters.
Despite successes in many use cases, previous approaches [9, 10, 11, 12] found
solutions lack of certifiable guarantees, which are unfortunately in need for
policy adoption in this safety-related application. Besides, their surrogate
approximations of connectivity may overlook important probabilistic scenarios.
This results in suboptimal quality of the generated plans. As application
domains for SMC solvers, this paper considers emergency shelter placement and
supply chain network management – two important stochastic connectivity
optimization problems.

It is challenging to solve SMC because of their highly intractable nature
($\text{NP}^{\text{PP}}$-complete) [13] – still intractable even with good
satisfiability solvers [14, 15, 16] and model counters [17, 18, 19, 20, 21,
22, 23]. Previous research on SMC solves either a special case or domain-
specific applications [24, 25, 26, 27, 28, 29, 30]. The special case is called
the Marginal Maximum-A-Posterior (MMAP) problem, whose decision version can be
formulated as a special case of SMC [31, 32, 33, 34, 35, 36] Both cases are
solved by optimizing the surrogate representations of the intractable model
counting in variational forms [37, 38], or via knowledge compilation [39, 36,
40] or via sample average approximation [41, 42, 43, 11, 44, 45, 46, 47].

Nevertheless, previous approaches either cannot quantify the quality of their
solutions, or offer one-sided guarantees, or offer guarantees which can be
arbitrarily loose. The lack of tight guarantees results in delayed policy
adoption in safety-related applications such as the stochastic connectivity
optimization considered in this paper. Second, optimizing surrogate objectives
without quantifying the quality of approximation leads to sub-optimal behavior
empirically. For example, previous stochastic connectivity optimization
solvers occasionally produce suboptimal plans because their surrogate
approximations overlook cases of significant probability. This problem is
amplified when combinatorial constraints are present.

We propose $\mathtt{XOR}$-$\mathtt{SMC}$,  a polynomial algorithm accessing
NP-oracles, to solve highly intractable SMC problems with constant
approximation guarantees. These guarantees hold with high (e.g. $>99\%$)
probability. The strong guarantees enable policy adoption in safety-related
domains and improve the empirical performance of SMC solving (e.g.,
eliminating sub-optimal behavior and providing constraint satisfaction
guarantees). The constant approximation means that the solver can correctly
decide the truth of an SMC formula if tightening or relaxing the bounds on the
model count by a multiplicative constant do not change its truth value. The
embedding algorithms allow us to find approximate solutions to beyond-$NP$ SMC
problems via querying $NP$ oracles. It expands the applicability of the state-
of-the-art SAT solvers to highly intractable problems.

The high-level idea behind $\mathtt{XOR}$-$\mathtt{SMC}$ is as follows.
Imagine a magic that randomly filters out half models (solutions) to an SAT
formula. Model counting can be approximated using this magic and an SAT solver
– we confirm the SAT formula has more than $2^{k}$ models if it is satisfiable
after applying this magic $k$ times. This magic can be implemented by
introducing randomized constraints. The idea is developed by researchers [48,
49, 50, 17, 18, 51, 52, 19, 20, 53]. In these works, model counting is
approximated with guarantees using polynomial algorithms accessing NP oracles.
$\mathtt{XOR}$-$\mathtt{SMC}$ notices such polynomial algorithms can be
encoded as SAT formulae. Hence, SAT-Modulo-Counting can be written as SAT-
Modulo-SAT (or equivalently SAT), when we embed the SAT formula compiled from
algorithms to solve model counting into SMC. The constant approximation
guarantee also carries.

We evaluate the performance of $\mathtt{XOR}$-$\mathtt{SMC}$ on real-world
stochastic connectivity optimization problems. In particular, we consider
applied problems of emergency shelter placement and supply chain management.
For the shelter placement problem, our $\mathtt{XOR}$-$\mathtt{SMC}$ finds
high-quality shelter assignments with less computation time and better quality
than competing baselines. For wheat supply chain management, the solutions
found by our $\mathtt{XOR}$-$\mathtt{SMC}$ are better than those found by
baselines and are close to the optimal solutions.

Our contributions can be summarized as follows:

  * •

We propose $\mathtt{XOR}$-$\mathtt{SMC}$ to solve highly intractable
Satisfiability Modulo Counting problems with a polynomial number of queries to
NP-oracles.

  * •

We prove the proposed $\mathtt{XOR}$-$\mathtt{SMC}$ method has a high
probability constant approximation guarantee.

  * •

In experiments, we consider two important SMC problems: emergency shelter
allocation and robust supply chain design. Our $\mathtt{XOR}$-$\mathtt{SMC}$
method requires less time to find better solutions compared to several
baselines.

##  2 Preliminaries

###  2.1 Satisfiability Modulo Theories

Satisfiability Modulo Theory (SMT) determines the SATisfiability (SAT) of a
Boolean formula, which contains predicates whose truth values are determined
by the background theory. SMT represents a line of successful efforts to build
general-purpose logic reasoning engines, encompassing complex expressions
containing bit vectors, real numbers, integers, and strings, etc [54]. Over
the years, many good SMT solvers are built, such as the Z3 [55, 56] and cvc5
[57]. They play a crucial role in automated theorem proving, program analysis
[58], program verification [59], and software testing [60].

###  2.2 Model Counting and Probabilistic Inference

Model counting computes the number of models (i.e., satisfying variable
assignments) to an SAT formula. Consider a Boolean formula $f(\mathbf{x})$,
where the input $\mathbf{x}$ is a vector of Boolean variables, and the output
$f$ is also Boolean. When we use 0 to represent false and 1 to represent true,
$\sum_{x}f(\mathbf{x})$ computes the model count. Model counting is closely
related to probabilistic inference and machine learning because the marginal
inference on a wide range of probabilistic models can be formulated as a
weighted model counting problem [61, 62].

Exact approaches for probabilistic inference and model counting are often
based on knowledge compilation [63, 21, 64, 65]. Approximate approaches
include Variational methods and sampling. Variational methods [27, 66, 67, 68,
69] use tractable forms to approximate a complex probability distribution. Due
to a tight relationship between counting and sampling [49], sampling-based
approaches are important for model counting. Importance sampling-based
techniques such as SampleSearch [70] is able to provide lower bounds. Markov
Chain Monte Carlo is asymptotically accurate. However, they cannot provide
guarantees except for a limited number of cases [71, 72]. The authors of [73,
74, 75] transform weighted integration into optimization queries using extreme
value distribution, which today is often called the “Gumbel trick” [76, 77].

###  2.3 XOR Counting

There is an interesting connection between model counting and solving
satisfiability problems subject to randomized XOR constraints. To illustrate
this, hold $\mathbf{x}$ at $\mathbf{x}_{0}$, suppose we would like to know if
$\sum_{\mathbf{y}\in\mathcal{Y}}f(\mathbf{x}_{0},\mathbf{y})$ exceeds $2^{q}$.
Consider the SAT formula:

$\displaystyle
f(\mathbf{x}_{0},\mathbf{y})\wedge\mathtt{XOR}_{1}(\mathbf{y})\wedge\ldots\wedge\mathtt{XOR}_{q}(\mathbf{y}).$
(1)

Here, $\mathtt{XOR}_{1},\ldots,\mathtt{XOR}_{q}$ are randomly sampled XOR
constraints. $\mathtt{XOR}_{i}(\mathbf{y})$ is the logical XOR or the parity
of a randomly-sampled subset of variables from $\mathbf{y}$. In other words,
$\mathtt{XOR}_{i}(\mathbf{y})$ is true if and only if an odd number of these
randomly sampled variables in the subset are true.

1 Randomly sample
$\mathtt{XOR}_{1}(\mathbf{y}),\ldots,\mathtt{XOR}_{q}(\mathbf{y})$;

2 if
_$f(\mathbf{x}_{0},\mathbf{y})\wedge\mathtt{XOR}_{1}(\mathbf{y})\wedge\ldots\wedge\mathtt{XOR}_{q}(\mathbf{y})$
is satisfiable_ then

3 return _True_ ;

4else

5 return _False_ ;

6 end if

Algorithm 1 $\mathtt{XOR}$-$\mathtt{Binary}$ ($f$, $\mathbf{x}_{0}$, $q$)

Formula (1) is likely to be satisfiable if more than $2^{q}$ different
$\mathbf{y}$ vectors render $f(\mathbf{x}_{0},\mathbf{y})$ true. Conversely,
Formula (1) is likely to be unsatisfiable if $f(\mathbf{x}_{0},\mathbf{y})$
has less than $2^{q}$ satisfying assignments. The significance of this fact is
that it essentially transforms model counting (beyond NP) into satisfiability
problems (within NP). An intuitive explanation of why this fact holds is that
each satisfying assignment $\mathbf{y}$ has 50% chance to satisfy a randomly
sampled XOR constraint. In other words, each XOR constraint “filters out” half
satisfying assignments. For example, the number of models satisfying
$f(\mathbf{x}_{0},\mathbf{y})\wedge\mathtt{XOR}_{1}(\mathbf{y})$ is
approximately half of that satisfying $f(\mathbf{x}_{0},\mathbf{y})$.
Continuing this chain of reasoning, if $f(\mathbf{x}_{0},\mathbf{y})$ has more
than $2^{q}$ solutions, there are still satisfying assignments left after
adding $q$ XOR constraints; hence formula (1) is likely satisfiable. The
reverse direction can be reasoned similarly. The precise mathematical argument
of the constant approximation is in Lemma 1.

######  Lemma 1.

[49, 17, 78] Given Boolean function $f(\mathbf{x}_{0},y)$ as defined above,

  * •

If $\sum_{\mathbf{y}}f(\mathbf{x}_{0},\mathbf{y})\geq 2^{q_{0}}$, then for any
$q\leq q_{0}$, with probability $1-\frac{2^{c}}{(2^{c}-1)^{2}}$,
$\mathtt{XOR}$-$\mathtt{Binary}$ $(f,\mathbf{x}_{0},q-c)$ returns true.

  * •

If $\sum_{\mathbf{y}}f(\mathbf{x}_{0},\mathbf{y})\leq 2^{q_{0}}$, then for any
$q\geq q_{0}$, with probability $1-\frac{2^{c}}{(2^{c}-1)^{2}}$,
$\mathtt{XOR}$-$\mathtt{Binary}$ $(w,\theta_{0},q+c)$ returns false.

This idea of transforming model counting problems into SAT problems subject to
randomized constraints is rooted in Leslie Valiant’s seminal work on unique
SAT [48, 49] and has been developed by a rich line of work [50, 17, 18, 51,
52, 19, 20, 53]. This idea has recently gathered momentum thanks to the rapid
progress in SAT solving [79, 16]. The contribution of this proposal extends
the success of SAT solvers to problems with even higher complexity, namely,
$\text{NP}^{\text{PP}}$-complete SMC problems.

##  3 Problem Formulation

Satisfiability Modulo Counting (SMC) is Satisfiability Modulo Theory (SMT)
[54] with model counting as the background theory. A canonical definition of
the SMC problem is to determine if there exists
$\mathbf{x}=(x_{1},\dots,x_{n})\in\mathcal{X}=\\{0,1\\}^{n}$ and
$\mathbf{b}=(b_{1},\dots,b_{k})\in\\{0,1\\}^{k}$ that satisfies the formula:

$\small\phi(\mathbf{x},\mathbf{b}),b_{i}\Leftrightarrow\left(\sum_{\mathbf{y}_{i}\in\mathcal{Y}_{i}}f_{i}(\mathbf{x},\mathbf{y}_{i})\geq
2^{q_{i}}\right),\forall i\in\\{1..,k\\}.$ (2)

Here each $b_{i}$ is a Boolean predicate that is true if and only if the
corresponding model count exceeds a threshold. Bold symbols (i.e.,
$\mathbf{x}$, $\mathbf{y}_{i}$ and $\mathbf{b}$) are vectors of Boolean
variables. $\phi,f_{1},\dots,f_{k}$ are Boolean functions (i.e., their input
is Boolean vectors, and their outputs are also Boolean). We use $0$ to
represent false and $1$ to represent true. Hence $\sum f_{i}$ computes the
number of satisfying assignments (model counts) of $f_{i}$. The directions of
the inequalities do not matter much because one can always negate each
$f_{i}$.

Our $\mathtt{XOR}$-$\mathtt{SMC}$ algorithm obtains the constant approximation
guarantee to the following slightly relaxed SMC problems. The problem
$\mathtt{SMC}(\phi,f_{1},\dots,f_{k},$ $q_{1},\dots,q_{k})$ finds a satisfying
assignment $(\mathbf{x},\mathbf{b})$ for:

$\displaystyle\phi(\mathbf{x},\mathbf{b})$
$\displaystyle\wedge\left[b_{1}\Rightarrow\left(\sum_{\mathbf{y}_{1}\in\mathcal{Y}_{1}}f_{1}(\mathbf{x},\mathbf{y}_{1})\geq
2^{q_{1}}\right)\right]\dots\wedge\left[b_{k}\Rightarrow\left(\sum_{\mathbf{y}_{k}\in\mathcal{Y}_{k}}f_{k}(\mathbf{x},\mathbf{y}_{k})\geq
2^{q_{k}}\right)\right].$ (3)

The only difference compared to the full-scale problem in Eq. (2)) is the
replacement of $\Leftrightarrow$ with $\Rightarrow$. This change allows us to
derive a concise constant approximation bound. We also mention that all the
applied SMC problems considered in this paper can be formulated in this
relaxed form.

##  4 The $\mathtt{XOR}$-$\mathtt{SMC}$ Algorithm

The key motivation behind our proposed $\mathtt{XOR}$-$\mathtt{SMC}$ algorithm
is to notice that Algorithm 1 itself can be written as a Boolean formula due
to the Cook-Levin reduction. When we embed this Boolean formula into Eq. (3),
the Satisfiability-Modulo-Counting problem translates into a Satisfiability-
Modulo-SAT problem, or equivalently, an SAT problem. This embedding also
ensures a constant approximation guarantee (see Theorem 2).

To illustrate the high-level idea, let us consider replacing each
$\sum_{\mathbf{y}_{i}\in\mathcal{Y}_{i}}f_{i}(\mathbf{x},\mathbf{y}_{i})\geq
2^{q_{i}}$ in Eq. (3) with formula

$\displaystyle
f_{i}(\mathbf{x},\mathbf{y}_{i})\wedge\mathtt{XOR}_{1}(\mathbf{y}_{i})\wedge\ldots\wedge\mathtt{XOR}_{q_{i}}(\mathbf{y}_{i}).$
(4)

We denote the previous equation (4) as
$\gamma(f_{i},\mathbf{x},q_{i},\mathbf{y}_{i})$. This replacement results in
the Boolean formula:

$\displaystyle\phi(\mathbf{x},\mathbf{b})\wedge$
$\displaystyle\left[b_{1}\Rightarrow\gamma(f_{1},\mathbf{x},q_{1},\mathbf{y}_{1})\right]\wedge\dots\wedge\left[b_{k}\Rightarrow\gamma(f_{k},\mathbf{x},q_{k},\mathbf{y}_{k})\right].$
(5)

Figure 1: $\mathtt{XOR}$-$\mathtt{SMC}$ (Algorithm 2) solves the intractable
model counting with satisfiability problems subject to randomized XOR
constraints and obtains constant approximation guarantees for SMC.

We argue that the satisfiability of formula (5) should be closely related to
that of formula (3) due to the connection between model counting and
satisfiability testing subject to randomized constraints (discussed in Section
2.3). To see this, Eq. (5) is satisfiable if and only if there exists
$(\mathbf{x},\mathbf{b},\mathbf{y}_{1},\ldots,\mathbf{y}_{k})$ that render Eq.
(5) true (notice $\mathbf{y}_{1},\ldots,\mathbf{y}_{k}$ are also its
variables). Suppose
$\texttt{SMC}(\phi,f_{1},\dots,f_{k},q_{1}+c,\dots,q_{k}+c)$ is satisfiable
(a.k.a., Eq. (3) is satisfiable when $q_{i}$ is replaced with $q_{i}+c$). Let
$(\mathbf{x},\mathbf{b})$ be a satisfying assignment. For any $b_{i}=1$ (true)
in $\mathbf{b}$, we must have
$\sum_{\mathbf{y_{i}}\in\mathcal{Y}_{i}}f_{i}(\mathbf{x},\mathbf{y}_{{i}})\geq
2^{q_{i}+c}$. This implies with a good chance, there exists a $\mathbf{y}_{i}$
that renders $\gamma(f_{i},\mathbf{x},q_{i},\mathbf{y}_{i})$ true. This is due
to the discussed connection between model counting and SAT solving subject to
randomized constraints. Hence
$b_{i}\Rightarrow\gamma(f_{i},\mathbf{x},q_{i},\mathbf{y}_{i})$ is true. For
any $b_{i}=0$ (false), the previous equation is true by default. Combining
these two facts and $\phi(\mathbf{x},\mathbf{b})$ is true, we see Eq. (5) is
true.

Conversely, suppose
$\texttt{SMC}(\phi,f_{1},\dots,f_{k},q_{1}-c,\dots,q_{k}-c)$ is not
satisfiable. This implies for every $(\mathbf{x},\mathbf{b})$, either
$\phi(\mathbf{x},\mathbf{b})$ is false, or there exists at least one $j$ such
that $b_{j}$ is true, but
$\sum_{\mathbf{y}_{{j}}\in\mathcal{Y}_{{j}}}f_{j}(\mathbf{x},\mathbf{y}_{{j}})<2^{q_{j}-c}$.
The first case implies Eq. (5) is false under the assignment. For the second
case,
$\sum_{\mathbf{y}_{{j}}\in\mathcal{Y}_{{j}}}f_{j}(\mathbf{x},\mathbf{y}_{{j}})<2^{q_{j}-c}$
implies with a good chance there is no $\mathbf{y}_{j}$ to make
$\gamma(f_{j},\mathbf{x},q_{j},\mathbf{y}_{j})$ true. Combining these two
facts, with a good chance Eq. (5) is not satisfiable.

In practice, to reduce the error probability the determination of the model
count needs to rely on the majority satisfiability status of a series of
equations (4) (instead of a single one). Hence we develop Algorithm 2, which
is a little bit more complex than the high-level idea discussed above. The
idea is still to transform the highly intractable SMC problem into solving an
SAT problem of its polynomial size, while ensuring a constant approximation
guarantee. Fig. 1 displays the encoding of Algorithm 2. We can see the core is
still to replace the intractable model count with satisfiability problems
subject to randomized constraints. We prove $\mathtt{XOR}$-$\mathtt{SMC}$ has
a constant approximation guarantee in Theorem 2. We leave the implementation
of $\mathtt{XOR}$-$\mathtt{SMC}$ in the Appendix B.

1 $T\leftarrow\lceil\frac{(n+k)\ln 2-\ln\eta}{\alpha(c,k)}\rceil$;

2

3for  _$t=1\text{ to }T$_ do

4 for  _$i=1\text{ to }k$_ do

5 $\psi_{i}^{(t)}\leftarrow f_{i}(\mathbf{x},\mathbf{y}_{i}^{(t)})$;

6 for  _$j=1\text{ to }q_{i}$_ do

7
$\psi_{i}^{(t)}\leftarrow\psi_{i}^{(t)}\wedge\mathtt{XOR}_{j}(\mathbf{y}_{i}^{(t)})$;

8

9 end for

10 $\psi_{i}^{(t)}\leftarrow\psi_{i}^{(t)}\vee\neg b_{i}$;

11

12 end for

13 $\psi_{t}\leftarrow\psi_{1}^{(t)}\wedge\dots\wedge\psi_{k}^{(t)}$;

14

15 end for

16$\phi^{*}\leftarrow\phi\wedge\mathtt{Majority}(\psi_{1},\ldots,\psi_{T})$ ;

17 if  _there exists
$(\mathbf{x},\mathbf{b},\\{\mathbf{y}_{i}^{(1)}\\}_{i=1}^{k},\ldots,\\{\mathbf{y}_{i}^{(T)}\\}_{i=1}^{k})$
that satisfies $\phi^{*}$_ then

18 return _True_ ;

19else

20 return _False_ ;

21 end if

Algorithm 2 $\mathtt{XOR}$-$\mathtt{SMC}$ ($\phi$, $\\{f_{i}\\}_{i=1}^{k}$,
$\\{q_{i}\\}_{i=1}^{k}$, $\eta$, $c$)

######  Theorem 2.

Let $0<\eta<1$ and $c\geq\log(k+1)+1$. Select $T=\lceil{((n+k)\ln
2-\ln\eta)}/{\alpha(c,k)}\rceil$, we have

  * •

Suppose there exists $\mathbf{x}_{0}\in\\{0,1\\}^{n}$ and
$\mathbf{b}_{0}\in\\{0,1\\}^{k}$, such that
$\mathtt{SMC}(\phi,f_{1},\dots,f_{k},q_{1}+c,\dots,q_{k}+c)$ is true. In other
words,

$\phi(\mathbf{x}_{0},\mathbf{b}_{0})\wedge\left(\bigwedge_{i=1}^{k}\left(b_{i}\Rightarrow\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x}_{0},\mathbf{y}_{i})\geq
2^{q_{i}+c}\right)\right),$

Then algorithm $\mathtt{XOR}$-$\mathtt{SMC}$ ($\phi$, $\\{f_{i}\\}_{i=1}^{k}$,
$\\{q_{i}\\}_{i=1}^{k}$, $T$) returns true with probability greater than
$1-\eta$.

  * •

Contrarily, suppose
$\mathtt{SMC}(\phi,f_{1},\dots,f_{k},q_{1}-c,\dots,q_{k}-c)$ is not
satisfiable. In other words, for all $\mathbf{x}\in\\{0,1\\}^{n}$ and
$\mathbf{b}\in\\{0,1\\}^{k}$,

$\displaystyle\neg\left(\phi(\mathbf{x},\mathbf{b})\wedge\left(\bigwedge_{i=1}^{k}\left(b_{i}\Rightarrow\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x},\mathbf{y}_{i})\geq
2^{q_{i}-c}\right)\right)\right),$

then $\mathtt{XOR}$-$\mathtt{SMC}$ ($\phi$, $\\{f_{i}\\}_{i=1}^{k}$,
$\\{q_{i}\\}_{i=1}^{k}$, $T$) returns false with probability greater than
$1-\eta$.

###### Proof.

Claim 1: Suppose there exists
$\mathbf{x}_{0}=[x_{1},\dots,x_{n}]\in\\{0,1\\}^{n}$ and
$\mathbf{b}_{0}=[b_{1},\dots,b_{k}]\in\\{0,1\\}^{k}$, such that

$\displaystyle\phi(\mathbf{x}_{0},\mathbf{b}_{0})\wedge\left(\bigwedge_{i=1}^{k}\left(b_{i}\Rightarrow\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x}_{0},\mathbf{y}_{i})\geq
2^{q_{i}+c}\right)\right)$ (6)

holds true. Denote $k_{0}$ as the number of non-zero bits in $\mathbf{b}_{0}$.
Without losing generality, suppose those non-zero bits are the first $k_{0}$
bits, i.e., $b_{1}=b_{2}=\dots=b_{k_{0}}=1$ and $b_{i}=0,\forall i>k_{0}$.
Then Eq. 6 can be simplified to:

$\displaystyle\phi(\mathbf{x}_{0},\mathbf{b}_{0})\wedge\left(\bigwedge_{i=1}^{k_{0}}\left(\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x}_{0},\mathbf{y}_{i})\geq
2^{q_{i}+c}\right)\right)$ (7)

Consider the Boolean formula $\psi_{t}$ defined in the
$\mathtt{XOR}$-$\mathtt{SMC}$ algorithm (choosing any $t\in\\{1,\dots,T\\}$).
$\psi_{t}$ can be simplified by substituting the values of $\mathbf{x}_{0}$
and $\mathbf{b}_{0}$. After simplification:

$\displaystyle\psi_{t}=$
$\displaystyle\left(f_{1}(\mathbf{x}_{0},\mathbf{y}_{1}^{(t)})\wedge\mathtt{XOR}_{1}(\mathbf{y}_{1}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{1}}(\mathbf{y}_{1}^{(t)})\right)\wedge\dots$
$\displaystyle\wedge\left(f_{k_{0}}(\mathbf{x}_{0},\mathbf{y}_{k_{0}}^{(t)})\wedge\mathtt{XOR}_{1}(\mathbf{y}_{k_{0}}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{k_{0}}}(\mathbf{y}_{k_{0}}^{(t)})\right)$

Denote
$\gamma_{i}=\left(f_{i}(\mathbf{x}_{0},\mathbf{y}_{i}^{(t)})\wedge\mathtt{XOR}_{1}(\mathbf{y}_{i}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{i}}(\mathbf{y}_{i}^{(t)})\right)$.
Observing that $\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x}_{0},\mathbf{y}_{i})\geq
2^{q_{i}+c},\forall i=1,\dots,k_{0}$, according to Lemma 1, with probability
at least $1-\frac{2^{c}}{(2^{c}-1)^{2}}$, there exists $\mathbf{y}_{i}^{(t)}$,
such that $(\mathbf{x}_{0},\mathbf{y}_{i}^{(t)})$ renders $\gamma_{i}$ true.
The probability that $\psi_{t}$ is true under
$(\mathbf{x}_{0},\mathbf{b}_{0},\mathbf{y}_{1}^{(t)},\dots,\mathbf{y}_{k}^{(t)})$
is:

$\displaystyle\mathbb{P}((\mathbf{x}_{0},\mathbf{b}_{0},\mathbf{y}_{1}^{(t)},\dots,\mathbf{y}_{k}^{(t)})\text{
renders }\psi_{t}\text{ true})$ $\displaystyle=~{}$
$\displaystyle\mathbb{P}\left(\bigwedge_{i=1}^{k_{0}}((\mathbf{x}_{0},\mathbf{y}_{i}^{(t)})\text{
renders }\gamma_{i}\text{ false})\right)$ $\displaystyle=~{}$ $\displaystyle
1-\mathbb{P}\left(\bigvee_{i=1}^{k_{0}}((\mathbf{x}_{0},\mathbf{y}_{i}^{(t)})\text{
renders }\gamma_{i}\text{ false})\right)$ $\displaystyle\geq~{}$
$\displaystyle
1-\sum_{i=1}^{k_{0}}\mathbb{P}((\mathbf{x}_{0},\mathbf{y}_{i}^{(t)})\text{
renders }\gamma_{i}\text{ false})$ $\displaystyle\geq~{}$ $\displaystyle
1-\frac{k_{0}2^{c}}{(2^{c}-1)^{2}}$ $\displaystyle\geq~{}$ $\displaystyle
1-\frac{k2^{c}}{(2^{c}-1)^{2}}.$

Define $\Gamma_{t}$ as a binary indicator variable where

$\displaystyle\Gamma_{t}=\begin{cases}1&\text{if
$(\mathbf{x}_{0},\mathbf{b}_{0},\mathbf{y}_{1}^{(t)},\dots,\mathbf{y}_{k}^{(t)})$
renders $\psi_{t}$ true}\\\ 0&\text{otherwise}\end{cases}$

Therefore $\mathbb{P}(\Gamma_{t}=0)\leq\frac{k2^{c}}{(2^{c}-1)^{2}}$.
$\mathbb{P}(\Gamma_{t}=0)<\frac{1}{2}$ when $c\geq\log_{2}(k+1)+1$.
$\mathtt{XOR}$-$\mathtt{SMC}$ returns true if the majority of
$\psi_{t},t=1,\dots,T$ are true; that is, $\sum_{t}\Gamma_{t}\geq\frac{T}{2}$.
Let’s define

$\displaystyle\alpha(c,k)$
$\displaystyle=D\left(\frac{1}{2}\|\frac{k2^{c}}{(2^{c}-1)^{2}}\right)=\frac{1}{2}\ln\frac{(2^{c}-1)^{2}}{k2^{c+1}}+\left(1-\frac{1}{2}\right)\ln\frac{2(2^{c}-1)^{2}}{(2^{c}-1)^{2}-k2^{c+1}}$

When $c\geq\log_{2}(k+1)+1$, observing that $\alpha(c,k)>0$ , we can apply the
Chernoff-Hoeffding theorem to obtain:

$\displaystyle\mathbb{P}\left(\sum_{t=1}^{T}\Gamma_{t}\geq\frac{T}{2}\right)$
$\displaystyle=1-\mathbb{P}\left(\sum_{t=1}^{T}\Gamma_{t}<\frac{T}{2}\right)\geq
1-e^{-\alpha(c,k)T}$

For $T\geq\lceil\frac{((n+k)\ln
2-\ln\eta)}{\alpha(c,k)}\rceil\geq\frac{-\ln\eta}{\alpha(c,k)}$, it follows
that $e^{-\alpha(c,k)T}\leq\eta$. Therefore, with a probability at least
$1-\eta$, we have $\sum_{t}\Gamma_{t}\geq\frac{T}{2}$. In this scenario,
$\mathtt{XOR}$-$\mathtt{SMC}$ ($\phi$, $\\{f_{i}\\}_{i=1}^{k}$,
$\\{q_{i}\\}_{i=1}^{k}$, $T$) returns true as it discovers $\mathbf{x}_{0}$,
$\mathbf{b}_{0}$, $(\mathbf{y}_{1}^{(t)},\dots,\mathbf{y}_{k}^{(t)})$, for
which the majority of Boolean formulae in $\\{\psi_{t}\\}_{t=1}^{T}$ are true.

Claim 2: Suppose for all $\mathbf{x}\in\\{0,1\\}^{n}$ and
$\mathbf{b}\in\\{0,1\\}^{k}$,

$\displaystyle\neg\left(\phi(\mathbf{x},\mathbf{b})\wedge\left(\bigwedge_{i=1}^{k}\left(b_{i}\Rightarrow\sum_{\mathbf{y}_{i}}f_{i}(\mathbf{x},\mathbf{y}_{i})\geq
2^{q_{i}-c}\right)\right)\right)$

Consider a fixed $\mathbf{x}_{1}$ and $\mathbf{b}_{1}$, we show the previous
condition with high probability renders most $\psi_{t}$ false in Algorithm 2.
We can prove that the probability is sufficiently low such that
$\mathtt{XOR}$-$\mathtt{SMC}$ will return false with a high probability after
examining all $\mathbf{x}$ and $\mathbf{b}$. The detailed proof is left to the
Appendix A. ∎

##  5 Experiment: Locate Emergency Shelters

In this section, we show our $\mathtt{XOR}$-$\mathtt{SMC}$ finds better
shelter location assignments (more emergency vacate paths) than baseline
approaches (in Fig. 3). It also needs less computational time (in Table 1).

###  5.1 Problem Formulation

Disasters such as hurricanes and floods continue to endanger millions of
lives. Shelters are safe zones that protect residents from possible damage,
and evacuation routes are the paths from resident zones toward shelter areas.
To enable the timely evacuation of resident zones, picking a set of shelter
locations with sufficient routing from resident areas should be considered.
Given the unpredictability and chaos during natural disasters, it is crucial
to guarantee multiple paths rather than a single path from residential areas
to shelters. This ensures that even if one route is obstructed, residents have
alternative paths to safety areas.

Current methods [80, 81] has considered finding those shelters locations that
have at least one single path from a residential area and model the problem of
finding an optimal assignment of shelters using Mixed Integer Programming
(MIP) framework [82, 83]. However, those proposed methods cannot be
generalized to solve the problem that requires sufficient alternative routes
from each residential area to shelters, primarily because counting the number
of paths is intractable. This complexity renders it nearly impossible to
address directly through MIP, especially on large-scale problems. Our proposed
$\mathtt{XOR}$-$\mathtt{SMC}$ can identify shelter location assignments that
ensure each residential area has a designated number of alternative paths with
high probability.

Figure 2: Example assignment of shelters that guarantee sufficient alternative
paths from the resident areas, at Hawaii Island. Every orange dot corresponds
to shelters and the green dot indicates a resident area.

Given a map $G=(V,E)$ where nodes in $V=\\{v_{1},\dots,v_{N}\\}$ represents
$N$ areas and $e=(v_{i},v_{j})\in E$ indicates a road from $v_{i}$ to $v_{j}$.
$N$ and $M$ denote the number of nodes and edges correspondingly. Given a
subset of nodes $R=\\{v_{r_{1}},\dots,v_{r_{k}}\\}\subseteq V$ indicates the
residential areas, the task is to choose at most $m$ nodes as shelters from
the rest nodes, such that the number of routes that can reach a shelter from
each residential area is maximized. Fig. 2 gives an example with $m=4$
shelters and there are many roads connecting the resident area to those
shelters.

SMC Formulation. We transform this optimization problem into a decision
problem by binary searching the thresholds on the number of paths. The
decision problem decides there are at least $2^{q_{r}}$ paths connecting any
residential area with a shelter. The assigned shelters is represented by a
vector $\mathbb{b}=(b_{1},\ldots,b_{n})\in\\{0,1\\}^{n}$, where $b_{i}=1$
implies node $v_{i}$ is chosen as shelter. Let
$\phi(\mathbf{b})=\left(\sum_{i=1}^{n}b_{i}\right)\leq m$ represent there are
at most $m$ shelters. Let $f(v_{r},v_{s},E^{\prime})$ be an indicator function
that returns one if and only if the selected edges $E^{\prime}$ form a path
from $v_{r}$ to $v_{s}$. The whole formula is:

$\displaystyle\phi(\mathbf{b}),b_{i}\Rightarrow$
$\displaystyle\left(\sum_{v_{s}\in S,E^{\prime}\subseteq
E}f(v_{r},v_{s},E^{\prime})\geq 2^{q_{r}}\right)\text{for }1\leq i\leq n.$

We leave the details implementation of $f(v_{r},v_{s},E^{\prime})$ in the
Appendix C.1.

###  5.2 Empirical Experiment Analysis

Experiment Setting. We crawl the real-world dataset from the Hawaii Statewide
GIS Program website. We extract the real Hawaii map with those major roads and
manually label those resident areas on the map. We create problems of
different scales by subtracting different sub-regions from the map. Fig. 2 is
an example of a labeled dataset with possible shelter assignments. 3 major
resident areas are picked as $R$, and set $m=5$.

Table 1: $\mathtt{XOR}$-$\mathtt{SMC}$ takes less empirical running time than baselines, for fining shelter location assignments over different graphs. Graph size is the number of nodes in the graph. | Graph Size
---|---
| $N=121$ | $N=183$ | $N=388$
$\mathtt{XOR}$-$\mathtt{SMC}$ (ours) | $\mathbf{0.04}h$ | $\mathbf{0.11}h$ | $\mathbf{0.16}h$
Gibbs-LS | $0.56h$ | $0.66h$ | $6.97h$
QuickSampler-LS | $0.31h$ | $0.29h$ | $0.62h$
Unigen-LS | $0.08h$ | $0.07h$ | $0.42h$

In terms of baselines, we consider the local search algorithm with shelter
locations as the state and the number of paths between shelters and resident
areas as the heuristic. Due to the intractability of path counting in our
formulation, the heuristic is approximated by querying sampling oracles. 1)
Gibbs sampling-based [84] Local Search (Gibbs-LS). 2) Uniform SAT sampler-
based [85] Local Search (Unigen-LS). 3) Quick Sampler-based [86] Local Search
(Quick-LS). Each baseline runs 5 times, and the best result is included.

For the evaluation metrics, 1) we consider the empirical running time when
reaching an assignment that guarantees each resident area has at least $2^{4}$
paths to a shelter, the result is shown in Table 1. 2) By gradually increasing
the path count threshold $q_{r}$, we can find the best shelter location that
guarantees the most number of paths. Given a time limit of 5 hours, we
evaluate different algorithms as in Fig. 3.

Result Analysis. In terms of running time analysis,
$\mathtt{XOR}$-$\mathtt{SMC}$ takes less empirical running time than baselines
for finding shelter location assignments over different graphs. Furthermore,
we evaluate the quality of the identified shelter locations. Our
$\mathtt{XOR}$-$\mathtt{SMC}$ finds shelter locations that have a higher
number of connecting paths than those found by baselines.

Figure 3: $\mathtt{XOR}$-$\mathtt{SMC}$ finds better shelter locations with
way more paths from the resident area to the chosen shelters than competing
baselines, across different graphs.

##  6 Experiment: Robust Supply Chain Design

In this section, we show our $\mathtt{XOR}$-$\mathtt{SMC}$ resolves robust
supply chain design problems against stochastic events. We show the solutions
found by $\mathtt{XOR}$-$\mathtt{SMC}$ outperform baselines and are near the
optimal ones in real-world applications.

###  6.1 Problem Formulation

Supply chain management founds its importance in logistics, operations
research, and economics. The essence of supply chain management is to
coordinate and integrate the flow of products, information, and finances
across companies to enhance overall performance and maximize value to the
consumer. Its importance is underscored by the increasing complexity of global
business environments, where even minor inefficiencies can result in
significant cost implications and jeopardize competitiveness.

Existing works in supply chain optimization often gravitate towards
mathematical programming approaches. However, when delving into more complex
scenarios involving stochastic events—such as uncertain demand or supply
disruptions—the task becomes considerably more intricate. Specifically,
formulating counting constraints, like guaranteeing a certain amount of
supplies across multiple vendors under stochastic events, is notably
intractable. These complexities necessitate innovative approaches that can
capture the randomness and dynamism inherent in real-world supply chain
systems without sacrificing optimality.

Suppose we have a supply chain of $N$ suppliers, they form a supply-demand
network $(V,E)$ where each node $v\in V$ represents a supplier and edges $e\in
E$ represent supply-to-demand trades. Each supplier $v$ plays two roles: as a
vendor to downstream suppliers, and as a buyer from upstream suppliers. There
are $M$ different kinds of products in the network. To guarantee substantial
production, supplier $v$ should order necessary raw materials from upstream
suppliers in advance. The trade between vendor $u$, and buyer $v$ costs
$f(u,v)$ for $v$. The capacity (amount of goods) of the trade between $u$ and
$v$ is $c(u,v)$. $v$ has total budget of $B(v)$ to get all his materials
ready. Due to unpredictable factors such as natural disasters, equipment
failure and etc., the trade between $u$ and $v$ may fail. Denote
$\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{L})\in\\{0,1\\}^{L}$ as the
state of $L$ different stochastic events, where $\theta_{l}=1$ indicates event
$l$ occurs. The task is to determine the optimal trading suppliers to maximize
expected production, accounting for stochastic influences.

SMC Formulation. The optimization problem can be formulated as a decision
problem similar to the previous section. The trading plan can be denoted as
$\\{x_{s,v}\\}_{s,v\in V}$, where $x_{s,v}=1$ indicates $v$ orders raw
material from $s$. For every node $v$, the decision problem decides whether

$\displaystyle\bigwedge_{v}\bigwedge_{d\in
D(v)}\mathbb{E}_{\boldsymbol{\theta}}\left(\sum_{s\in
S(d)}c(s,v)I(x_{s,v},\boldsymbol{\theta})\right)\geq 2^{q_{v,d}}$ (8)

where $D(v)$ denotes the set of materials required by $v$, $S(d)$ represents
the set of suppliers that produce raw material $d$, $2^{q_{v,d}}$ is a
threshold measuring the amount of product $d$ can be acquired by $v$, and
indicator function $I(x_{s,v},\boldsymbol{\theta})$ is evaluated to one if the
random events do not affect the trade between $s$ and $v$. Eq. (8) also needs
to be accompanied by a budget constraint per supplier. We leave the details of
implementation in Appendix C.2.

###  6.2 Empirical Experiment Analysis

Experiment Setting. We evaluate our algorithm on the wheat supply chain
network from [87]. A detailed introduction is shown in the appendix. Although
the cost of trade, the capacity of transportation, raw material demand and
etc. can be directly extracted, it lacks a model of stochastic events due to
the limitation of their algorithm. We generated 10 stochastic events (see
Appendix C.2) overall supply-demand edges that jointly determine the trading,
which makes the expectation intractable by enumeration.

For the baseline, we used Gibbs-CD [88] and Belief Propagation with
Contrastive Divergence (BP-CD) [89] as two alternative methods for estimating
the maximum expected productivity for every node which corresponds to their
trading strategy. As for the metric, the expected total productivity (tons of
bread) of all markets is evaluated in Table 2. For a better comparison, two
synthetic networks with sizes of 20 and 100, in which the parameters are
generated with a Gaussian distribution fitted from the real-world data.

Result Analysis. The productivity of bread is used in the evaluation.
$\mathtt{XOR}$-$\mathtt{SMC}$ finds the most optimized supply chain network
among all methods. It also produces solutions close to the best possible ones.

Table 2:  Expectation of the total bread production (unit: ton) from the supply network. $\mathtt{XOR}$-$\mathtt{SMC}$ proposes the best bread production supply chain that achieves the maximum total productivity and is close to the best possible solutions. | Supply Network Size
---|---
|  $20$ (Synth) |  $44$ (Real) |  $100$ (Synth)
$\mathtt{XOR}$-$\mathtt{SMC}$ (ours) | $\mathbf{947.61}$ | $\mathbf{3467.82}$ | $\mathbf{4753.78}$
Gibbs-CD | $911.43$ | $3217.70$ | $4190.39$
BP-CD | $917.66$ | $3325.51$ | $4374.66$
Best Possible | $1055.22$ | $4223.23$ | $6046.48$

##  7 Conclusion

We presented $\mathtt{XOR}$-$\mathtt{SMC}$, an algorithm with polynomial
approximation guarantees to solve the highly intractable Satisfiability Modulo
Counting (SMC) problems. SMC fuses statistical inference and symbolic
reasoning. Solving it presents unique challenges. Prior work on SMC solving
offers no or loose guarantees and may find suboptimal solutions.
$\mathtt{XOR}$-$\mathtt{SMC}$ transforms the intractable SMC problem into
satisfiability problems by replacing intricate model counting with SAT
formulae subject to randomized XOR constraints. $\mathtt{XOR}$-$\mathtt{SMC}$
also obtains constant approximation guarantees on the solutions obtained. SMC
solvers offer useful tools for many real-world problems at the nexus of
symbolic and statistical AI. Extensive experiments on SMC problems for social
good demonstrate that $\mathtt{XOR}$-$\mathtt{SMC}$ outperforms other
approaches in solution quality and also runs efficiently.

## Acknowledgements

This research was supported by NSF grants IIS-1850243, CCF-1918327.

## References

  * Freuder and O’Sullivan [2023] Eugene C. Freuder and Barry O’Sullivan, editors.  _AAAI-23 Constraint Programming and Machine Learning Bridge Program_ , 2023.
  * Belle et al. [2023] Vaishak Belle, Michael Fisher, Xiaowei Huang, Masoumeh Mansouri, Albert Merono-Penuela, Sriraam Natarajan, and Efi Tsamoura, editors.  _Neuro-symbolic AI for Agent and Multi-Agent Systems (NeSyMAS) Workshop at AAMAS-23_ , 2023.
  * Munawar et al. [2023] Asim Munawar, Jon Lenchner, Francesca Rossi, Lior Horesh, Alexander Gray, and Murray Campbell, editors.  _IBM Neuro-Symbolic AI Workshop 2023 – Unifying Statistical and Symbolic AI_ , 2023.
  * Wu et al. [2015] Xiaojian Wu, Daniel Sheldon, and Shlomo Zilberstein.  Fast combinatorial algorithm for optimizing the spread of cascades.  In _Proceedings of the 24th International Joint Conference on Artificial Intelligence (IJCAI)_ , pages 2655–2661, 2015.
  * Dilkina and Gomes [2010] Bistra Dilkina and Carla P. Gomes.  Solving connected subgraph problems in wildlife conservation.  In _CPAIOR_ , pages 102–116, 2010.
  * Israeli and Wood [2002] Eitan Israeli and R Kevin Wood.  Shortest-path network interdiction.  _Networks: An International Journal_ , 40(2):97–111, 2002.
  * Kempe et al. [2005] David Kempe, Jon Kleinberg, and Éva Tardos.  Influential nodes in a diffusion model for social networks.  In _Automata, languages and programming_ , pages 1127–1138. Springer, 2005.
  * Almeida et al. [2019] R. Almeida, Qinru Shi, Jonathan M. Gomes-Selman, Xiaojian Wu, Yexiang Xue, H. Angarita, N. Barros, B. Forsberg, R. García-Villacorta, S. Hamilton, J. Melack, M. Montoya, Guillaume Perez, S. Sethi, C. Gomes, and A. Flecker.  Reducing greenhouse gas emissions of amazon hydropower with strategic dam planning.  _Nature Communications_ , 10, 2019.
  * Williams and Snyder [2005] Justin C. Williams and Stephanie A. Snyder.  Restoring habitat corridors in fragmented landscapes using optimization and percolation models.  _Environmental Modeling and Assessment_ , 10(3):239–250, 2005.
  * Conrad et al. [2012] John Conrad, Carla P. Gomes, Willem-Jan van Hoeve, Ashish Sabharwal, and Jordan F. Suter.  Wildlife corridors as a connected subgraph problem.  _Journal of Environmental Economics and Management_ , 63(1), 2012.
  * Sheldon et al. [2010a] Daniel Sheldon, Bistra Dilkina, Adam Elmachtoub, Ryan Finseth, Ashish Sabharwal, Jon Conrad, Carla Gomes, David Shmoys, William Allen, Ole Amundsen, and William Vaughan.  Maximizing the spread of cascades using network design.  In _Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence_ , pages 517–526, 2010a.
  * Wu et al. [2014] Xiaojian Wu, Daniel Sheldon, and Shlomo Zilberstein.  Stochastic network design in bidirected trees.  In _Advances in Neural Information Processing Systems_ , pages 882–890, 2014.
  * Park and Darwiche [2004] James D. Park and Adnan Darwiche.  Complexity results and approximation strategies for map explanations.  _J. Artif. Int. Res._ , 2004.
  * Biere et al. [2009] Armin Biere, Marijn J. H. Heule, Hans van Maaren, and Toby Walsh, editors.  _Handbook of Satisfiability_ , volume 185 of _Frontiers in Artificial Intelligence and Applications_.  IOS Press, February 2009.  ISBN 978-1-58603-929-5.
  * Rossi et al. [2006] Francesca Rossi, Peter van Beek, and Toby Walsh.  _Handbook of Constraint Programming (Foundations of Artificial Intelligence)_.  Elsevier Science Inc., New York, NY, USA, 2006.  ISBN 0444527265.
  * Braunstein et al. [2005] A. Braunstein, M. Mézard, and R. Zecchina.  Survey propagation: an algorithm for satisfiability.  _Random Struct. Algorithms_ , 27:201–226, 2005.
  * Gomes et al. [2006] Carla P. Gomes, Ashish Sabharwal, and Bart Selman.  Model counting: A new strategy for obtaining good bounds.  In _Proceedings of the 21st National Conference on Artificial Intelligence_ , 2006.
  * Ermon et al. [2013a] Stefano Ermon, Carla P. Gomes, Ashish Sabharwal, and Bart Selman.  Taming the curse of dimensionality: Discrete integration by hashing and optimization.  In _Proceedings of the 30th International Conference on Machine Learning, ICML_ , 2013a.
  * Achlioptas and Theodoropoulos [2017] Dimitris Achlioptas and Panos Theodoropoulos.  Probabilistic model counting with short xors.  In _International Conference on Theory and Applications of Satisfiability Testing_ , pages 3–19. Springer, 2017.
  * Chakraborty et al. [2013] Supratik Chakraborty, Kuldeep S. Meel, and Moshe Y. Vardi.  A scalable and nearly uniform generator of sat witnesses.  In _Proceedings of the 25th International Conference on Computer Aided Verification_ , 2013.
  * Kisa et al. [2014] Doga Kisa, Guy Van den Broeck, Arthur Choi, and Adnan Darwiche.  Probabilistic sentential decision diagrams.  In _Proceedings of the 14th International Conference on Principles of Knowledge Representation and Reasoning (KR)_ , 2014.
  * Cheng et al. [2012] Qiang Cheng, Feng Chen, Jianwu Dong, Wenli Xu, and Alexander Ihler.  Approximating the sum operation for marginal-map inference.  In _Conf. on Artificial Intelligence (AAAI)_ , pages 1882–1887. July 2012.
  * Gogate and Dechter [2012] Vibhav Gogate and Rina Dechter.  Importance sampling-based estimation over and/or search spaces for graphical models.  _Artificial Intelligence_ , 184-185:38 – 77, 2012.
  * Belanger and McCallum [2016] David Belanger and Andrew McCallum.  Structured prediction energy networks.  In _International Conference on Machine Learning_ , pages 983–992, 2016.
  * Welling and Teh [2003] Max Welling and Yee Whye Teh.  Approximate inference in boltzmann machines.  _Artificial Intelligence_ , 143(1):19–50, 2003\.
  * Yedidia et al. [2001] Jonathan S Yedidia, William T Freeman, and Yair Weiss.  Generalized belief propagation.  In _Advances in neural information processing systems_ , pages 689–695, 2001.
  * Wainwright and Jordan [2008] Martin J. Wainwright and Michael I. Jordan.  Graphical models, exponential families, and variational inference.  _Found. Trends Mach. Learn._ , 1(1-2):1–305, January 2008.
  * Fang et al. [2015] Fei Fang, Peter Stone, and Milind Tambe.  When security games go green: Designing defender strategies to prevent poaching and illegal fishing.  In _IJCAI_ , 2015.
  * Conitzer and Sandholm [2006] Vincent Conitzer and Tuomas Sandholm.  Computing the optimal strategy to commit to.  In _Proc. 7th ACM Conference on Electronic Commerce (EC)_ , pages 82–90, 2006.
  * Sheldon et al. [2010b] Daniel Sheldon, Bistra N. Dilkina, Adam N. Elmachtoub, Ryan Finseth, Ashish Sabharwal, Jon Conrad, Carla P. Gomes, David B. Shmoys, William Allen, Ole Amundsen, and William Vaughan.  Maximizing the spread of cascades using network design.  In _UAI_ , 2010b.
  * Marinescu et al. [2014] Radu Marinescu, Rina Dechter, and Alexander T. Ihler.  AND/OR search for marginal MAP.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI_ , 2014.
  * Liu and Ihler [2013] Qiang Liu and Alexander T. Ihler.  Variational algorithms for marginal MAP.  _Journal of Machine Learning Research_ , 14, 2013.
  * Mauá and de Campos [2012] Denis Deratani Mauá and Cassio Polpo de Campos.  Anytime marginal MAP inference.  In _Proceedings of the 29th International Conference on Machine Learning, ICML_ , 2012.
  * Jiang et al. [2011] Jiarong Jiang, Piyush Rai, and Hal Daumé III.  Message-passing for approximate MAP inference with latent variables.  In _Advances in Neural Information Processing Systems 24_ , 2011.
  * Lee et al. [2016] Junkyu Lee, Radu Marinescu, Rina Dechter, and Alexander T. Ihler.  From exact to anytime solutions for marginal MAP.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, AAAI_ , 2016.
  * Ping et al. [2015] Wei Ping, Qiang Liu, and Alexander T. Ihler.  Decomposition bounds for marginal MAP.  In _Advances in Neural Information Processing Systems 28_ , 2015.
  * Liu and Ihler [2012] Qiang Liu and Alexander Ihler.  Belief propagation for structured decision making.  In _Uncertainty in Artificial Intelligence (UAI)_ , pages 523–532. AUAI Press, Corvallis, Oregon, August 2012.
  * Kiselev and Poupart [2014] Igor Kiselev and Pascal Poupart.  Policy optimization by marginal-map probabilistic inference in generative models.  In _Proceedings of the 2014 international conference on Autonomous agents and multi-agent systems_ , pages 1611–1612. Citeseer, 2014.
  * Choi et al. [2022] YooJung Choi, Tal Friedman, and Guy Van den Broeck.  Solving marginal map exactly by probabilistic circuit transformations.  In _Proceedings of the 25th International Conference on Artificial Intelligence and Statistics (AISTATS)_ , mar 2022.
  * Mei et al. [2018] Jun Mei, Yong Jiang, and Kewei Tu.  Maximum A posteriori inference in sum-product networks.  In _AAAI_ , pages 1923–1930. AAAI Press, 2018.
  * Kleywegt et al. [2002] Anton J. Kleywegt, Alexander Shapiro, and Tito Homem-de Mello.  The sample average approximation method for stochastic discrete optimization.  _SIAM Journal on Optimization_ , 12:479–502, February 2002\.
  * Shapiro [2003] Alexander Shapiro.  Monte carlo sampling methods.  _Handbooks in operations research and management science_ , 10:353–425, 2003.
  * Swamy and Shmoys [2006] Chaitanya Swamy and David B. Shmoys.  Approximation algorithms for 2-stage stochastic optimization problems.  In _FSTTCS 2006: Foundations of Software Technology and Theoretical Computer Science_ , pages 5–19, 2006.
  * Dyer and Stougie [2006] Martin Dyer and Leen Stougie.  Computational complexity of stochastic programming problems.  _Mathematical Programming_ , 106(3):423–432, 2006\.
  * Wu et al. [2017] Xiaojian Wu, Akshat Kumar, Daniel Sheldon, and Shlomo Zilberstein.  Robust optimization for tree-structured stochastic network design.  In _Proceedings of the 31st Conference on Artificial Intelligence_ , 2017.
  * Xue et al. [2015] Shan Xue, Alan Fern, and Daniel Sheldon.  Scheduling conservation designs for maximum flexibility via network cascade optimization.  _J. Artif. Intell. Res. (JAIR)_ , 2015.
  * Verweij et al. [2003] Bram Verweij, Shabbir Ahmed, Anton J. Kleywegt, George Nemhauser, and Alexander Shapiro.  The sample average approximation method applied to stochastic routing problems: A computational study.  _Computational Optimization and Applications_ , 24(2-3):289–333, 2003.
  * Valiant and Vazirani [1986] L. G. Valiant and V. V. Vazirani.  Np is as easy as detecting unique solutions.  _Theoretical Computer Science_ , 47:85–93, 1986.
  * Jerrum et al. [1986] M.R. Jerrum, L.G. Valiant, and V.V. Vazirani.  Random generation of combinatorial structures from a uniform distribution.  _Theoretical Computer Science_ , 43:169–188, 1986.  ISSN 0304-3975.
  * Gomes et al. [2007] Carla P Gomes, Ashish Sabharwal, and Bart Selman.  Near-uniform sampling of combinatorial spaces using xor constraints.  In _Advances in Neural Information Processing Systems_ , 2007.
  * Ermon et al. [2013b] Stefano Ermon, Carla P. Gomes, Ashish Sabharwal, and Bart Selman.  Embed and project: Discrete sampling with universal hashing.  In _Advances in Neural Information Processing Systems (NIPS)_ , 2013b.
  * Kuck et al. [2019] Jonathan Kuck, Tri Dao, Shengjia Zhao, Burak Bartan, Ashish Sabharwal, and Stefano Ermon.  Adaptive hashing for model counting.  In _Conference on Uncertainty in Artificial Intelligence_ , 2019.
  * Chakraborty et al. [2014] Supratik Chakraborty, Daniel J. Fremont, Kuldeep S. Meel, Sanjit A. Seshia, and Moshe Y. Vardi.  Distribution-aware sampling and weighted model counting for sat.  In _AAAI_ , 2014.
  * Barrett et al. [2009] Clark Barrett, Roberto Sebastiani, Sanjit Seshia, and Cesare Tinelli.  _Satisfiability Modulo Theories_ , chapter 26, pages 825–885.  2009\.
  * De Moura and Bjørner [2008] Leonardo De Moura and Nikolaj Bjørner.  Z3: An efficient smt solver.  In _International conference on Tools and Algorithms for the Construction and Analysis of Systems_ , pages 337–340. Springer, 2008.
  * Bjørner et al. [2019] Nikolaj Bjørner, Leonardo de Moura, Lev Nachmanson, and Christoph M Wintersteiger.  Programming z3.  _Engineering Trustworthy Software Systems: 4th International School, SETSS 2018, Chongqing, China, April 7–12, 2018, Tutorial Lectures 4_ , pages 148–201, 2019.
  * Barbosa et al. [2022] Haniel Barbosa, Clark Barrett, Martin Brain, Gereon Kremer, Hanna Lachnitt, Makai Mann, Abdalrhman Mohamed, Mudathir Mohamed, Aina Niemetz, Andres Nötzli, et al.  cvc5: A versatile and industrial-strength smt solver.  In _International Conference on Tools and Algorithms for the Construction and Analysis of Systems_ , pages 415–442. Springer, 2022.
  * Feser et al. [2020] John Feser, Sam Madden, Nan Tang, and Armando Solar-Lezama.  Deductive optimization of relational data storage.  _Proceedings of the ACM on Programming Languages_ , 4(OOPSLA):1–30, 2020.
  * VK et al. [2022] Hari Govind VK, Sharon Shoham, and Arie Gurfinkel.  Solving constrained horn clauses modulo algebraic data types and recursive functions.  _Proc. ACM Program. Lang._ , 6(POPL):1–29, 2022\.
  * De Moura and Bjørner [2007] Leonardo De Moura and Nikolaj Bjørner.  Efficient e-matching for smt solvers.  In _Automated Deduction–CADE-21: 21st International Conference on Automated Deduction Bremen, Germany, July 17-20, 2007 Proceedings 21_ , pages 183–198. Springer, 2007.
  * Chavira and Darwiche [2008] Mark Chavira and Adnan Darwiche.  On probabilistic inference by weighted model counting.  _Artificial Intelligence_ , 172(6-7):772–799, 2008.
  * Xue et al. [2016a] Yexiang Xue, Zhiyuan Li, Stefano Ermon, Carla P. Gomes, and Bart Selman.  Solving marginal map problems with np oracles and parity constraints.  In _Proceedings of the 29th Annual Conference on Neural Information Processing Systems (NIPS)_ , 2016a.
  * Darwiche and Marquis [2002] Adnan Darwiche and Pierre Marquis.  A knowledge compilation map.  _J. Artif. Int. Res._ , 2002.
  * Choi et al. [2013] Arthur Choi, Doga Kisa, and Adnan Darwiche.  Compiling probabilistic graphical models using sentential decision diagrams.  In _Proceedings of the 12th European Conference on Symbolic and Quantitative Approaches to Reasoning with Uncertainty (ECSQARU)_ , 2013.
  * Xue et al. [2012] Yexiang Xue, Arthur Choi, and Adnan Darwiche.  Basing decisions on sentences in decision diagrams.  In _Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence (AAAI)_ , 2012.
  * Wainwright et al. [2003] Martin J. Wainwright, Tommi S. Jaakkola, and Alan S. Willsky.  Tree-reweighted belief propagation algorithms and approximate ML estimation by pseudo-moment matching.  In _Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics_ , 2003.
  * Sontag et al. [2008] David Sontag, Talya Meltzer, Amir Globerson, Tommi Jaakkola, and Yair Weiss.  Tightening lp relaxations for map using message passing.  In _UAI_ , pages 503–510, 2008.
  * Hazan and Shashua [2010] T. Hazan and A. Shashua.  Norm-product belief propagation: Primal-dual message-passing for approximate inference.  _Information Theory, IEEE Transactions on_ , 56(12):6294–6316, Dec 2010.  ISSN 0018-9448.  doi: 10.1109/TIT.2010.2079014.
  * Flerova et al. [2011] Natalia Flerova, Er Ihler, Rina Dechter, and Lars Otten.  Mini-bucket elimination with moment matching.  In _In NIPS Workshop DISCML_ , 2011.
  * Gogate and Dechter [2007] Vibhav Gogate and Rina Dechter.  Samplesearch: A scheme that searches for consistent samples.  In _Artificial Intelligence and Statistics_ , pages 147–154. PMLR, 2007.
  * Jerrum and Sinclair [1996] Mark Jerrum and Alistair Sinclair.  The markov chain monte carlo method: an approach to approximate counting and integration.  _Approximation algorithms for NP-hard problems_ , pages 482–520, 1996\.
  * Madras [2002] N.N. Madras.  _Lectures on Monte Carlo Methods_.  American Mathematical Society, 2002.  ISBN 0821829785.
  * Papandreou and Yuille [2010] George Papandreou and Alan L Yuille.  Gaussian sampling by local perturbations.  In _Advances in Neural Information Processing Systems 23_ , pages 1858–1866. 2010.
  * Hazan and Jaakkola [2012] Tamir Hazan and Tommi S. Jaakkola.  On the partition function and random maximum a-posteriori perturbations.  In _Proceedings of the 29th International Conference on Machine Learning, ICML_ , 2012.
  * Balog et al. [2017] Matej Balog, Nilesh Tripuraneni, Zoubin Ghahramani, and Adrian Weller.  Lost relatives of the gumbel trick.  In _Proceedings of the 34th International Conference on Machine Learning, ICML_ , pages 371–379, 2017.
  * Papandreou and Yuille [2011] George Papandreou and Alan L Yuille.  Perturb-and-map random fields: Using discrete optimization to learn and sample from energy models.  In _2011 International Conference on Computer Vision_ , pages 193–200. IEEE, 2011.
  * Jang et al. [2017] Eric Jang, Shixiang Gu, and Ben Poole.  Categorical reparameterization with gumbel-softmax.  In _ICLR (Poster)_. OpenReview.net, 2017.
  * Ermon et al. [2013c] Stefano Ermon, Carla Gomes, Ashish Sabharwal, and Bart Selman.  Taming the curse of dimensionality: Discrete integration by hashing and optimization.  In _International Conference on Machine Learning_ , pages 334–342. PMLR, 2013c.
  * Maneva et al. [2005] E. Maneva, Elchanan Mossel, and M. Wainwright.  A new look at survey propagation and its generalizations.  _ArXiv_ , cs.CC/0409012, 2005.
  * Bayram and Yaman [2018] Vedat Bayram and Hande Yaman.  Shelter location and evacuation route assignment under uncertainty: A benders decomposition approach.  _Transportation science_ , 52(2):416–436, 2018\.
  * Amideo et al. [2019] A Esposito Amideo, Maria Paola Scaparra, and Kathy Kotiadis.  Optimising shelter location and evacuation routing operations: The critical issues.  _European Journal of Operational Research_ , 279(2):279–295, 2019.
  * Liang et al. [2019] Bian Liang, Dapeng Yang, Xinghong Qin, and Teresa Tinta.  A risk-averse shelter location and evacuation routing assignment problem in an uncertain environment.  _International journal of environmental research and public health_ , 16(20):4007, 2019.
  * Hansuwa et al. [2022] Sweety Hansuwa, Usha Mohan, and Viswanath Kumar Ganesan.  Scenario-based stochastic shelter location-allocation problem with vulnerabilities for disaster relief network design.  _European Journal of Industrial Engineering_ , 16(5):507–530, 2022.
  * Geman and Geman [1984] Stuart Geman and Donald Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 6(6):721–741, 1984.
  * Soos et al. [2020] Mate Soos, Stephan Gocht, and Kuldeep S. Meel.  Tinted, detached, and lazy cnf-xor solving and its applications to counting and sampling.  In _Proceedings of International Conference on Computer-Aided Verification (CAV)_ , 7 2020.
  * Dutra et al. [2018] Rafael Dutra, Kevin Laeufer, Jonathan Bachrach, and Koushik Sen.  Efficient sampling of sat solutions for testing.  In _Proceedings of the 40th international conference on software engineering_ , pages 549–559, 2018.
  * Zokaee et al. [2017] Shiva Zokaee, Armin Jabbarzadeh, Behnam Fahimnia, and Seyed Jafar Sadjadi.  Robust supply chain network design: an optimization model with real world application.  _Annals of Operations Research_ , 257:15–44, 2017.
  * Carreira-Perpinan and Hinton [2005] Miguel A Carreira-Perpinan and Geoffrey E Hinton.  On contrastive divergence learning.  In _Aistats_ , volume 10, pages 33–40. Citeseer, 2005.
  * Fan and Yexiang [2020] Ding Fan and Xue Yexiang.  Contrastive divergence learning with chained belief propagation.  In _International Conference on Probabilistic Graphical Models_ , pages 161–172. PMLR, 2020.
  * Jerrum et al. [2004] Mark Jerrum, Alistair Sinclair, and Eric Vigoda.  A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries.  _Journal of the ACM (JACM)_ , 51(4):671–697, 2004\.
  * Xue et al. [2016b] Yexiang Xue, Zhiyuan Li, Stefano Ermon, Carla P Gomes, and Bart Selman.  Solving marginal map problems with np oracles and parity constraints.  _Advances in Neural Information Processing Systems_ , 29, 2016b.

##  Appendix A Proofs of Theorem 2

###  A.1 Proof of Claim 2

Suppose for all $\mathbf{x}\in\\{0,1\\}^{n}$ and $\mathbf{b}\in\\{0,1\\}^{k}$,
there is

$\displaystyle\neg\left(\phi(\mathbf{x},\mathbf{b})\wedge\left(\bigwedge_{i=1}^{k}\left(b_{i}\Rightarrow\sum_{y_{i}}f_{i}(\mathbf{x},y_{i})\geq
2^{q_{i}-c}\right)\right)\right)$
$\displaystyle\Rightarrow\neg\phi(\mathbf{x},\mathbf{b})\vee\left(\bigvee_{i=1}^{k}\left(b_{i}\wedge\sum_{y_{i}}f_{i}(\mathbf{x},y_{i})\leq
2^{q_{i}-c}\right)\right).$

Consider a fixed $\mathbf{x}_{1}$ and $\mathbf{b}_{1}$, we will examine the
probability of not returning true when discovering
$(\mathbf{x}_{1},\mathbf{b}_{1})$.

If $\neg\phi(\mathbf{x}_{1},\mathbf{b}_{1})$ holds true, then the probability
of $\mathtt{XOR}$-$\mathtt{SMC}$ ($\phi$, $\\{f_{i}\\}_{i=1}^{k}$,
$\\{q_{i}\\}_{i=1}^{k}$, $T$) not returning true for discovering
$(\mathbf{x}_{1},\mathbf{b}_{1})$ is one.

Otherwise, when
$\bigvee_{i=1}^{k}\left(b_{i}\wedge\sum_{y_{i}}f_{i}(\mathbf{x}_{1},y_{i})\leq
2^{q_{i}-c}\right)$ is true, i.e,

$\displaystyle\exists
i,~{}s.t.~{}~{}b_{i}\wedge\sum_{y_{i}}f_{i}(\mathbf{x}_{1},y_{i})\leq
2^{q_{i}-c}.$ (9)

Let’s select any $t=1,\dots,T$. The Boolean formula $\psi_{t}$ can be
simplified by substituting the values of $\mathbf{x}_{1}$ and
$\mathbf{b}_{1}$:

$\displaystyle\psi_{t}=$
$\displaystyle\left(f_{1}(\mathbf{x}_{1},y_{1}^{(t)})\wedge\mathtt{XOR}_{1}(y_{1}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{1}}(y_{1}^{(t)})\right)\wedge\dots$
$\displaystyle\wedge\left(f_{k_{1}}(\mathbf{x}_{1},y_{k_{1}}^{(t)})\wedge\mathtt{XOR}_{1}(y_{k_{1}}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{k_{1}}}(y_{k_{1}}^{(t)})\right)$

where $k_{1}$ is the number of digits that are non-zero in $\mathbf{b}_{1}$,
and $1\leq k_{1}\leq k$.

According to Eq. (9), we define

$\gamma_{i}=\left(f_{i}(\mathbf{x}_{1},y_{i}^{(t)})\wedge\mathtt{XOR}_{1}(y_{i}^{(t)})\wedge\dots\wedge\mathtt{XOR}_{q_{i}}(y_{i}^{(t)})\right)$

where $\sum_{y_{i}^{(t)}}f_{i}(\mathbf{x}_{1},y_{i}^{(t)})\leq 2^{q_{i}-c}$.
Then according to Lemma 1, with probability at least
$1-\frac{2^{c}}{(2^{c}-1)^{2}}$, $(\mathbf{x}_{1},y_{i}^{(t)})$ renders
$\gamma_{i}$ false for all $y_{i}^{(t)}$.

Then the probability of $\psi_{t}$ being false under all
$(\mathbf{x}_{1},\mathbf{b}_{1},y_{1}^{(t)},\dots,y_{k}^{(t)})$ is

$\displaystyle\mathbb{P}((\mathbf{x}_{1},\mathbf{b}_{1},y_{1}^{(t)},\dots,y_{k}^{(t)})\text{
renders $\psi_{t}$ false for all $\\{y_{i}^{(t)}\\}_{i=1}^{k}$})$
$\displaystyle=~{}$
$\displaystyle\mathbb{P}\left(\bigvee_{i=1}^{k_{1}}(\mathbf{x}_{1},y_{i}^{(t)})\text{
renders $\gamma_{i}$ false for all $y_{i}^{(t)}$}\right)$
$\displaystyle\geq~{}$
$\displaystyle\mathbb{P}((\mathbf{x}_{1},y_{i}^{(t)})\text{ renders
$\gamma_{i}$ false for all $y_{i}^{(t)}$})$ $\displaystyle\geq~{}$
$\displaystyle 1-\frac{2^{c}}{(2^{c}-1)^{2}}$

Define $\Gamma_{t}$ as a binary indicator variable where

$\displaystyle\Gamma_{t}=\begin{cases}0&\text{if all
$(\mathbf{x}_{1},\mathbf{b}_{1},y_{1}^{(t)},\dots,y_{k}^{(t)})$ renders
$\psi_{t}$ false}\\\ 1&\text{otherwise}\end{cases}$

Upon discovering $\mathbf{x}_{1},\mathbf{b}_{1}$,
$\mathtt{XOR}$-$\mathtt{SMC}$ won’t return true if the majority of $\psi_{t}$,
$t=1,\dots,T$ are false under
$(\mathbf{x}_{1},\mathbf{b}_{1},y_{1}^{(t)},\dots,y_{k}^{(t)})$ for all
$\\{y_{i}^{(t)}\\}_{i=1}^{k}$, that is,
$\sum_{t=1}^{T}\Gamma_{t}\leq\frac{T}{2}$. By Chernoff–Hoeffding theorem,

$\displaystyle\mathbb{P}\left(\sum_{t=1}^{T}\Gamma_{t}\leq\frac{T}{2}\right)$
$\displaystyle=1-\mathbb{P}\left(\sum_{t=1}^{T}\Gamma_{t}>\frac{T}{2}\right)\geq
1-e^{-\alpha(c,1)T}$

For $T\geq\frac{(n+k)\ln 2-\ln\eta}{\alpha(c,1)}$, it follows that
$e^{-\alpha(c,1)T}\leq\frac{\eta}{2^{n+k}}$, In this case, with a probability
at least $1-\frac{\eta}{2^{n+k}}$, $\mathtt{XOR}$-$\mathtt{SMC}$ won’t return
true upon discovering $\mathbf{x}_{1}$ and $\mathbf{b}_{1}$.

Noted that $\mathtt{XOR}$-$\mathtt{SMC}$ returns false only if it won’t return
true upon discovering any $(\mathbf{x}_{1},\mathbf{b}_{1})$. Therefore, we can
write

$\displaystyle\mathbb{P}(\text{$\mathtt{XOR}$-$\mathtt{SMC}$ returns false})$
$\displaystyle=~{}$
$\displaystyle\mathbb{P}(\forall\mathbf{x},\mathbf{b},\text{$\mathtt{XOR}$-$\mathtt{SMC}$
won't return true upon $(\mathbf{x},\mathbf{b})$})$ $\displaystyle=~{}$
$\displaystyle
1-\mathbb{P}(\exists\mathbf{x}_{1},\mathbf{b}_{1},\text{$\mathtt{XOR}$-$\mathtt{SMC}$
returns true upon $(\mathbf{x}_{1},\mathbf{b}_{1})$})$ $\displaystyle\geq~{}$
$\displaystyle
1-\sum_{\mathbf{x}_{1},\mathbf{b}_{1}}\mathbb{P}(\text{$\mathtt{XOR}$-$\mathtt{SMC}$
returns true upon $(\mathbf{x}_{1},\mathbf{b}_{1})$})$ $\displaystyle\geq~{}$
$\displaystyle 1-\frac{\eta}{2^{n+k}}2^{n+k}$ $\displaystyle=~{}$
$\displaystyle 1-\eta.$

###  A.2 Proof of Theorem 2

Combining both Claim 1 and 2, we can choose the value $T=\frac{(n+k)\ln
2-\ln\eta}{\alpha(c,k)}$ such that
$T\geq\max\\{\frac{-\ln\eta}{\alpha(c,k)},\frac{(n+k)\ln
2-\ln\eta}{\alpha(c,1)}\\}$. This completes the proof.

##  Appendix B Implementation of $\mathtt{XOR}$-$\mathtt{SMC}$

Our code is attached to the supplementary material. It contains 1) the
implementation of our $\mathtt{XOR}$-$\mathtt{SMC}$ method 2) the list of
datasets, and 3) the implementation of several baseline algorithms.

We use CPLEX solver111https://www.ibm.com/products/ilog-cplex-optimization-
studio/cplex-optimizer version 22.1.1 as the NP oracle. We use Python package
boolexpr222http://www.boolexpr.org/index.html to process those parts involve
logic operations “$\lor,\land,\neg$”.

###  B.1 Computational Pipeline of $\mathtt{XOR}$-$\mathtt{SMC}$

#### Preprocessing

The preprocessing step converts the problem into a standard form as specified
in Eq. (3):

$\displaystyle\phi(\mathbf{x},\mathbf{b})$
$\displaystyle\wedge\left[b_{1}\Rightarrow\left(\sum_{\mathbf{y}_{1}\in\mathcal{Y}_{1}}f_{1}(\mathbf{x},\mathbf{y}_{1})\geq
2^{q_{1}}\right)\right]\dots\wedge\left[b_{k}\Rightarrow\left(\sum_{\mathbf{y}_{k}\in\mathcal{Y}_{k}}f_{k}(\mathbf{x},\mathbf{y}_{k})\geq
2^{q_{k}}\right)\right].$

The crucial aspects are to guarantee

  * •

$\forall i$, $\mathcal{Y}_{i}$ should be properly encoded such that
$\mathcal{Y}_{i}\subseteq\\{0,1\\}^{d_{i}}$ for some $d_{i}$. (See example in
Shelter Location Assignment.)

  * •

$\forall i$, $f_{i}:\mathcal{X}\times\mathcal{Y}_{i}\rightarrow\\{0,1\\}$ must
be a binary-output function. (See example in Robust Supply Chain Design.)

#### Formulate $f$ into a Boolean Formula

A binary-output function $f_{i}$ can be formulated into a Boolean formula
where $\mathbf{x}$ and $\mathbf{y}_{i}$ are seen as boolean variables and
$f_{i}(\mathbf{x},\mathbf{y}_{i})=1$ iff. the boolean formula is satisfiable
at $(\mathbf{x},\mathbf{y}_{i})$.

Although it seems straightforward to use SAT solvers, we found it more
efficient to transform the SAT into an equivalent Constraint Programming (CP)
problem and query the CPLEX solver after comparison. Due to the
interchangeability of CP and SAT, we will keep using the SAT formulation for
the following discussion.

#### Add XOR Constraints

Then we need to add random XOR constraints to the Boolean formula. As an
example, an XOR constraint for $\mathbf{y}\in\\{0,1\\}^{d}$ is

$\displaystyle\mathtt{XOR}(\mathbf{y})=y_{2}\oplus y_{3}\oplus y_{d}\oplus 1$

In this $\mathtt{XOR}(\mathbf{y})$, each literal $y_{j}$, $j=1,\dots,d$, and
the extra $1$ in the end (adding an extra $1$ can be seen as negating all
literals) has $1/2$ chance of appearance. Following this manner, we will
sample multiple XOR constraints for each $\mathbf{y}_{i}$. The number of XOR
constraints for $\mathbf{y}_{i}$ is determined by the logarithm of the
corresponding threshold $\log_{2}(2^{q_{i}})=q_{i}$.

After sampling all XOR constraints, we will have a boolean formula:

$\displaystyle(b_{1}\Rightarrow
f_{1}(\mathbf{x},\mathbf{y}_{1})\wedge\mathtt{XOR}_{1}(\mathbf{y}_{1})\wedge\ldots\wedge\mathtt{XOR}_{q_{1}}(\mathbf{y}_{1}))\wedge$
$\displaystyle(b_{2}\Rightarrow
f_{2}(\mathbf{x},\mathbf{y}_{2})\wedge\mathtt{XOR}_{2}(\mathbf{y}_{2})\wedge\ldots\wedge\mathtt{XOR}_{q_{2}}(\mathbf{y}_{2}))\wedge$
(10) $\displaystyle\cdots$

#### Repeat and Examine the Majority

Since all XOR constraints are sampled uniformly, the chance of Eq. (B.1) being
satisfiable when Eq. (3) holds true may not be significant. Therefore we can
enhance the probability by repeating experiments $T$ times and checking
whether the majority is satisfiable.

##  Appendix C Experiment Settings

All experiments run on a high-end server with two Sky Lake CPUs @ 2.60GHz, 24
Core, 96 GB RAM.

###  C.1 Shelter Location Assignment

As discussed in the main text, the essential decision problem decides there
are at least $2^{q_{r}}$ paths connecting any residential area with a shelter.
The assigned shelters is represented by a vector
$\mathbb{b}=(b_{1},\ldots,b_{n})\in\\{0,1\\}^{n}$, where $b_{i}=1$ implies
node $v_{i}$ is chosen as shelter. Let
$\phi(\mathbf{b})=\left(\sum_{i=1}^{n}b_{i}\right)\leq m$ represent there are
at most $m$ shelters. Let $I(v_{r},v_{s},E^{\prime})$ be an indicator function
that returns one if and only if the selected edges $E^{\prime}$ form a path
from $v_{r}$ to $v_{s}$. The whole formula is:

$\displaystyle\phi(\mathbf{b}),b_{i}\Rightarrow$
$\displaystyle\left(\sum_{v_{s}\in S,E^{\prime}\subseteq
E}I(v_{r},v_{s},E^{\prime})\geq 2^{q_{r}}\right)\text{for }1\leq i\leq n.$

#### Path Encoding

The path indicator function $I(\cdot)$ is implemented as follows. Any path
between 2 nodes can be mapped to a unique flow where the flow along the path
is 1 and others are 0. And one such unit integer flow can be represented by

$\mathbf{f}\in\\{(f_{e_{1}},f_{e_{2}},\dots,f_{e_{M}})|f_{e_{i}}\in\\{0,1\\},\text{
for }1\leq i\leq M\\},$

where $f_{e_{i}}=1$ if and only if the amount of flow on edge $e_{i}$ is 1.
Then we can encode every path on the map using $\mathbf{f}$. Then we can use a
flow indicator function instead in the implementation, where

$\displaystyle I(v_{r},v_{s},\mathbf{f})=\begin{cases}1&\text{$\mathbf{f}$
represents a path from $v_{r}$ to shelter $v_{s}$},\\\
0&\text{otherwise}.\end{cases}$

Then $I(\cdot)$ can be implemented by encoding flow constraints.

#### Shelter Location Encoding

The set of nodes with shelters $S$ can be further encoded by a binary vector

$\mathbf{x}_{S}=(x_{1},\dots,x_{N})\in\\{0,1\\}^{N},$

where $x_{i}=1$ represents a shelter built on node $v_{i}$, i.e., $v_{i}\in
S$. For example,

$\mathbf{x}_{S}=(1,1,1,\underbrace{0,\dots,0}_{\text{all $0$}})$ (11)

implies we place only 3 shelters on node $v_{1}$, $v_{2}$, and $v_{3}$.

Therefore, the shelter location problem can be formulated into an SMC:

$\displaystyle\left(\sum_{i=1}^{N}x_{i}\leq m\right)\wedge\bigwedge_{v_{r}\in
R}\left(\sum_{v_{s}\in S}\sum_{\mathbf{f}}I(v_{r},v_{s},\mathbf{f})\geq
2^{q_{r}}\right)$ (12)

If problem (12) holds true for some $S$, then the shelter assignment plan
meets the safety requirement. Thus in $\mathtt{XOR}$-$\mathtt{SMC}$, we are
actually solving an SAT problem:

$\displaystyle\left(\sum_{i=1}^{N}x_{i}\leq
m\right)\wedge\bigwedge_{v_{r}}\underbrace{\left(I(v_{r},v_{s},\mathbf{f})\wedge\mathtt{XOR}_{1}(\mathbf{f},\mathbf{x}_{S})\wedge\dots\wedge\mathtt{XOR}_{q_{r}}(\mathbf{f},\mathbf{x}_{S})\right)}_{\text{repeat
$T$ times and check the satisfiability of the majority}}$

where we determine whether there is such a $\mathbf{f}$ and $\mathbf{x}_{S}$
renders the boolean formula above true.

#### Dataset

The dataset is from Shelters at
Hawaii333https://iu.maps.arcgis.com/apps/Styler/index.html?appid=477fae419a014a3ba7c95cdcb39fd946.
The case study for this assignment is from Hawaii and involves Mauna Loa
volcano. This map offers an analysis of the most practical locations for
emergency shelters for displaced people within 15 minutes of their original
locations. Relationships between population, emergency shelters, and high-risk
lava flow zones on the islands are analyzed.

According to different choices of sub-regions, we extract 3 maps with 121,
186, and 388 nodes representing 3 difficulty levels. Those maps are included
in the code repository.

#### Evaluation Metric

Suppose we find a shelter assignment $S=\\{v_{s_{1}},\dots,v_{s_{m}}\\}$ for
residential areas $R=\\{v_{r_{1}},\dots,v_{r_{p}}\\}$.

The quality of the assignment is evaluated by the total paths from the
resident areas to the chosen shelters, which is computed as:

$\mathtt{\\#Path}(R,S)=\sum_{v_{r}\in R}\sum_{v_{s}\in S}\text{Number of paths
from $v_{r}$ to $v_{s}$}$

#### Choice of Baseline

For baseline methods, shelter locations are found by local search algorithms
defined as Alg. 3, where we search for $m$ best shelter locations. Function
randomNode($V$) represents initializing the starting search location randomly
from nodes $V$, and getNeighbors($v$) enumerates all neighbor nodes of $v$.
The heuristic function $h$ is defined as

$\displaystyle h(R,S)=~{}\sum_{v_{r}\in R}\sum_{v_{s}\in S}\text{Approx.
number of paths from $v_{r}$ to $v_{s}$}$

Due to the high computational cost of path counting, we instead use the
approximated number of paths as the heuristic suggesting the quality of a
shelter location. The equivalence of approximated counting and sampling has
been widely studied [90], so we implemented the approximated counter by
querying different sampling oracles including Gibbs Sampler, Unigen
Sampler444https://github.com/meelgroup/unigen, Quick
sampler555https://github.com/RafaelTupynamba/quicksampler, etc.

1 $S\leftarrow\\{\\}$;

2 for  _$i=1\textbf{ to }m$_ do

3 $v_{curr}\leftarrow$ randomNode($V$);

4 $score_{curr}\leftarrow h(R,\\{v_{curr}\\})$;

5 while  _True_ do

6 $neighbors\leftarrow$ getNeighbors($v_{curr}$);

7 $hasImproved\leftarrow False$;

8 for  _$v_{neigh}\textbf{ in }neighbors$_ do

9 $score_{neigh}\leftarrow h(R,\\{v_{neigh}\\})$;

10 if  _$score_{neigh} >score_{curr}$_ then

11 $v_{curr}\leftarrow v_{neigh}$;

12 $score_{curr}\leftarrow score_{neigh}$;

13 $hasImproved\leftarrow True$;

14 break;

15

16 end if

17

18 end for

19 if  _not $hasImproved$_ then

20 $S\leftarrow S\cup\\{v_{curr}\\}$;

21 break;

22

23 end if

24

25 end while

26

27 end for

return _$S$_ ;

Algorithm 3 Local-Search($G=(V,E)$, $R=\\{v_{r_{1}},\dots,v_{r_{p}}\\}$, $m$)

#### Hyper Parameters

For all experiments, residential nodes are picked at node $v_{r_{1}}=0$,
$v_{r_{2}}=10$, $v_{r_{3}}=20$. The maximum number of shelters is $m=5$.
Parameter $T$ (for repeat experiments) is set to be 2, which empirically works
well.

Table 1 shows the run time of finding a shelter assignment with a fixed
threshold, all results are generated when setting
$2^{q_{r_{1}}}=2^{q_{r_{2}}}=2^{q_{r_{3}}}=2^{4}$.

Due to the sensitivity of local search to the search state initialization. We
initialize the starting point randomly and repeat 5 times. Only those best are
picked for comparison.

###  C.2 Robust Supply Chain Design

As discussed in the main text Eq. (8), the core decision problem is defined as

$\displaystyle\bigwedge_{v}\bigwedge_{d\in
D(v)}\mathbb{E}_{\boldsymbol{\theta}}\left(\sum_{s\in
S(d)}c(s,v)I(x_{s,v},\boldsymbol{\theta})\right)\geq 2^{q_{v,d}}$

where $\\{x_{s,v}\\}_{s,v\in V}$ is the trading plan, $D(v)$ denotes the set
of materials required by $v$, $S(d)$ represents the set of suppliers that
produce raw material $d$, $2^{q_{v,d}}$ is a threshold measuring the amount of
product $d$ can be acquired by $v$, and indicator function
$I(x_{s,v},\boldsymbol{\theta})$ is evaluated to one if the random events do
not affect the trade between $s$ and $v$.

$\displaystyle I(x_{s,v},\boldsymbol{\theta})=\begin{cases}1&\text{$x_{s,v}=1$
and $\boldsymbol{\theta}$ doesn't break the trade,}\\\
0&\text{otherwise}.\end{cases}$

Besides, it also needs to be accompanied by a budget constraint per supplier.

$\displaystyle\bigwedge_{v}\left(\sum_{s}f(s,v)x_{s,v}\leq B(v)\right)$ (13)

Integrating (13) results in the SMC formulation givens a final form in
Equation 3.

#### Discretization

The expectation in Eq. (8) can be expanded as

$\displaystyle\mathbb{E}_{\boldsymbol{\theta}}\left(\sum_{s\in
S(d)}c(s,v)I(x_{s,v},\boldsymbol{\theta})\right)=~{}\sum_{\boldsymbol{\theta}}\sum_{s\in
S(d)}\mathbb{P}(\boldsymbol{\theta})c(s,v)I(x_{s,v},\boldsymbol{\theta})$

However, as required by Theorem 2, the counted term should be a binary output
function, so we need to discretize $\mathbb{P}(\boldsymbol{\theta})c(s,v)$ by
introducing a discretization scheme borrowed from [91].

Generally, for any weighted function
$w(x):\\{0,1\\}^{n}\rightarrow\mathbb{R}^{+}$, define the following embedding

$\displaystyle S(w,l)=\left\\{(x,y)|\forall 1\leq i\leq
l,\frac{w(x)}{M}\leq\frac{2^{i-1}}{2^{l}}\Rightarrow y_{i}=0\right\\}$

where $M=\max_{x}w(x)$.

######  Lemma 3 ([91]).

Let $w_{l}^{\prime}(x,y)$ be an indicator variable which is 1 if and only if
$(x,y)$ is in $S(w,l)$, i.e., $w_{l}^{\prime}(x,y)=\mathbf{1}_{(x,y)\in
S(w,l)}$, we claim that,

$\displaystyle\sum_{x}w(x)\leq\frac{M}{2^{l}}\sum_{(x,y)}w_{l}^{\prime}(x,y)\leq
2\sum_{x}w(x)+M2^{n-l}$

Lemma 3 provides a discretization method that can transform a general weighted
function into a sum of binary function. Applying this discretization method to
both $\mathbb{P}(\boldsymbol{\theta})$ and $c(s,v)$ results in a standard SMC
form as required.

#### Dataset

We evaluate our algorithm on the wheat supply chain network from [87], the
structure of the supply chain is shown in Fig. 4. For a better comparison, two
synthetic networks with sizes of 20 and 100, in which the parameters are
generated with a Gaussian distribution fitted from the real-world data. The
detailed Statistics of synthetic datasets for the supply network is collected
in Table 3.

Table 3: Statistics of synthetic datasets for the supply network. | #Wheat suppliers | #Flour factories | #Bread factories | #Markets | Total nodes in graph
---|---|---|---|---|---
20 (Synth) | 5 | 5 | 5 | 5 | 20
44 (Real) | 9 | 7 | 9 | 19 | 44
100 (Synth) | 25 | 25 | 25 | 25 | 100
Figure 4: Example of the wheat-to-bread supply chain network proposed by [87].
The example includes 3 wheat suppliers, 2 flour factories, 2 bread factories,
and 2 markets. A black solid arrow from $i$ to $j$ indicates a trade between
them ($x_{i,j}=1$ in the trading plan), and a dashed arrow from $p$ to $q$
represents an unchosen trade ($x_{p,q}=0$). Each edge also incorporates a cost
$f$ and capacity $c$ which are not shown in the figure. The red cross
indicates the edge is destroyed due to stochastic events (with a probability
specified by $\mathbb{P}(\boldsymbol{\theta})$).

#### Evaluation Metric

The trading plan is evaluated by the total expected "productivity" of markets,
which is equivalent to measuring the expected amount of raw material flows to
all markets. Formally defined by,

$\displaystyle\mathbb{E}_{\boldsymbol{\theta}}\left(\sum_{v_{m}}\sum_{v_{b}}c(v_{b},v_{m})I(x_{v_{b},v_{m}},\boldsymbol{\theta})\right)=~{}\sum_{\boldsymbol{\theta}}\sum_{v_{m}}\sum_{v_{b}}\mathbb{P}(\boldsymbol{\theta})c(v_{b},v_{m})I(x_{v_{b},v_{m}},\boldsymbol{\theta})$

where $v_{m}$ represents a market node and $v_{b}$ represents a bread factory
node.

#### Choice of Baseline

Baseline methods estimate the maximum expectation using Contrastive Divergence
Learning methods from samples. As shown in Table 2, we used Gibbs-CD [88] and
Belief Propagation with Contrastive Divergence (BP-CD) [89] as two alternative
methods for estimating the maximum expected productivity.

#### Hyper Parameters

Parameter $T$ (for repeat experiments) is set to be 2, which is sufficient to
show superior performance over baselines.

Both probability $P(\boldsymbol{\theta})$ and capacity functions $c(\cdot)$
are discretized into 4 bits, i.e., set $l=4$ in Lemma 3.

We consider 10 stochastic events on all trading edges in the network and all
events are assumed independent, i.e.,
$\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{10})$ and
$\mathbb{P}(\boldsymbol{\theta})=\mathbb{P}(\theta_{1})\mathbb{P}(\theta_{2})\cdots\mathbb{P}(\theta_{10})$.

The detailed network parameter configurations including stochastic event
definition are specified in the code repository.