# Minimizing the Outage Probability in a Markov Decision Process

Vincent Corlay and Jean-Christophe Sibel
Mitsubishi Electric R&D Centre Europe, Rennes, France. E-mail: {v.corlay,
<EMAIL_ADDRESS>
###### Abstract

Standard Markov decision process (MDP) and reinforcement learning algorithms
optimize the policy with respect to the expected gain. We propose an algorithm
which enables to optimize an alternative objective: the probability that the
gain is greater than a given value. The algorithm can be seen as an extension
of the value iteration algorithm. We also show how the proposed algorithm
could be generalized to use neural networks, similarly to the deep $Q$
learning extension of $Q$ learning.

###### Index Terms:

Markov decision process, value iteration, neural network, outage probability.

##  I Introduction

We consider an agent trying to learn the optimal sequence of actions to be
taken in an environment in order to maximize the total amount of reward. This
total amount of reward is called the gain. Standard approaches, both in the
framework of MDP and reinforcement learning, focus on the expected (average)
gain. Consequently, most algorithms in the literature are designed to maximize
this expected gain. For instance, all algorithms in the reference books [3]
and [5] optimize the expected gain. See Section I-A for connections with
recent works. As an illustration in the scope of a resource allocation
problem, a scheduler (the agent in this case) should decide which entity gets
the resource at each time step, see [4] for an example. In the standard
approach, the guideline is to minimize the average number of system failures.

However, one may be interested in ensuring that the gain is almost always
greater than a given value, i.e., optimizing the outage probability. This is
relevant for safety applications where some situations should be avoided with
the highest possible probability, even if the average performance is reduced.
In the case of a resource allocation problem, it can be relevant to minimize
the probability that the number of system failures is greater than a given
quantity. In information theory, the outage probability of a communication
channel, defined as the probability that a given information rate is not
supported, is an important metric.

Consequently, we propose an algorithm, inspired from existing MDP algorithms,
which enables to optimize the outage probability. We also explain how this
algorithm can be generalized to allow the use of neural networks, instead of
storing many values in a table.

###  I-A Connections with recent relevant works

A relevant fundamental work from the literature is [1]. In a nutshell, they
propose to work directly with random variables, and thus distributions, rather
than the expected gain. As a result, they consider the distributional Bellman
equation as an alternative to the standard Bellman equation:

$\displaystyle Z(s,a)=R(s,a)+\lambda Z(s^{\prime},a^{\prime}),$ (1)

where $Z$ is the gain random variable (denoted by $G$ in this paper) given
action $a$ in state $s$. This leads to a Monte Carlo algorithm (Algorithm 1 in
the paper) where atoms and probabilities are tracked in the iteration loop.

However, the objective in [1] is not the outage probability but still the
expected gain. There is also a key implementation difference: We focus on the
MDP approach, meaning that we assume the knowledge of the model (i.e., the
transition probabilities between the states). In [1], the practical part
(Section 4) considers the Monte Carlo approach without the model.
Consequently, the way we update the probabilities is different. We use e.g.,
(14) whereas it is based on Monte Carlo realizations $(s,a,r,s^{\prime})$ in
[1].

##  II Standard maximum expected gain approach

We consider an infinite-horizon MDP with a discounted sum-reward criterion as
presented in [3, Chap. 6], defined by states, actions, transition
probabilities, and rewards. The sum-reward criterion, called the gain, is
computed as

$\displaystyle G_{t}=\sum_{i=t}^{\infty}\lambda^{i-t}r_{i},$ (2)

where $r_{i}$ is the reward received at each time step $i$ and
$0\leq\lambda<1$ the discount.

In general, the objective in this framework is to find a policy $\pi$, i.e.,
the decision rule to be used at each time step, that maximizes the expected
gain $\mathbb{E}[G_{t}]$. Here, the expectation is applied with respect to the
randomness in the environment, modelled by transition probabilities between
the system states.

Accordingly, the value of a state $s\in S$ obtained with a given policy $\pi$,
where $S$ is the set of possible system states, is defined as the expected
gain given that the system is in the state $s$ at time $t$:

$\displaystyle v^{\pi}(s)=\mathbb{E}[G_{t}|s].$ (3)

Then, a policy $\pi^{*}$ is optimal if

$\displaystyle v^{\pi^{*}}(s)\geq v^{\pi}(s),\ \forall s\in S\text{ and
}\forall\pi.$ (4)

Under the optimal policy, the value of a state $s$ can be expressed via
Bellman equation as

$\displaystyle v^{\pi^{*}}(s)=\max_{a\in A_{s}}Q(s,a),$ (5)

where $A_{s}$ is the set of allowable actions in state $s$ and $Q(s,a)$,
called the $Q$-value, is defined as

$\displaystyle Q(s,a)=\mathbb{E}[G_{t}|s,a]=R(s,a)+\lambda\sum_{s_{j}\in
S}p(s_{j}|s,a)v^{\pi^{*}}(s_{j}),$ (6)

with $p(s_{j}|s,a)$ denoting the transition probability from state $s$ to
state $s_{j}$ given action $a$ and the average short-term reward $R(s,a)$ is
$R(s,a)=\sum_{j}p(s_{j}|s,a)r(s_{j},s)$, where $r(s_{j},s)$ is the reward
obtained when going from state $s$ to state $s_{j}$. Given the $Q$-values, the
optimal action in a state $s$ is obtained as

$\displaystyle a^{*}=\underset{a\in A_{s}}{\text{ arg max }}Q(s,a).$ (7)

A standard approach to compute the values $v^{\pi^{*}}(s)$ is to use the value
iteration algorithm [5, Chap 4.4]. It consists in computing $v_{n+1}(s)$ at
iteration $n+1$, $\forall s\in S$, as

$\displaystyle v_{n+1}(s)=\max_{a\in A_{s}}\ \\{R(s,a)+\lambda\sum_{s_{j}\in
S}p(s_{j}|s,a)v_{n}(s_{j})\\},$ (8)

where $v_{n}(s)$ converges to $v^{\pi^{*}}(s)$ as it is the unique fixed point
(see [3] for a proof). This step can be seen as a training phase.

Having the values of $v^{\pi^{*}}(s)$ $\forall s\in S$, the action to be taken
in a given state $s$ according to $\pi^{*}$ can be implemented by computing
(6) $\forall a\in A_{s}$ and then using (7). We refer to the obtained policy
as the maximum expected gain policy.

##  III Alternative proposed approach

###  III-A Problem statement

Instead of maximizing the expected gain, we propose to maximize the
probability that $G_{t}$ is greater than a given value $\alpha$. Since
$p(G_{t}>\alpha)=1-p(G_{t}\leq\alpha)$, this is equivalent to minimizing
$p(G_{t}<\alpha)$, commonly called outage probability. This outage probability
represents the risk to have a bad system outcome which is essential for safety
applications. Therefore, the value of a state $s$ with parameter $\alpha$ is
now defined as follows:

$\displaystyle v^{\pi}(s,\alpha)=p(G_{t}>\alpha|s).$ (9)

Accordingly, we also define the value of a state with parameter $\alpha$ given
an action $a\in A_{s}$:

$\displaystyle Q(s,a,\alpha)=p(G_{t}>\alpha|s,a).$ (10)

The new goal is thus to find a policy $\pi^{*}$ that maximizes the probability
that the gain $G_{t}$ is above an arbitrary value $\alpha$, i.e., such that:

$\displaystyle v^{\pi^{*}}(s,\alpha)\geq v^{\pi}(s,\alpha),\ \forall s\in
S\text{ and }\forall\pi.$ (11)

This latter policy $\pi^{*}$ is called the alternative policy.

###  III-B Example

We consider the recycling robot problem, presented in [5, Chap. 3], as a toy
example to illustrate the difference between the two objectives. The problem
is the following: A mobile robot running on a battery should collect empty
cans. It can be in two battery states
$S=\\{s_{1}=\text{``low"},s_{2}=\text{``high"}\\}$. For both states, the
possible actions are $A=\\{a_{1}=\text{``search", }a_{2}=\text{``wait"
},a_{3}=\text{``recharge" }\\}$. With the action “wait”, the robot stays in
the same state with probability $p(s_{i}|s_{i},a_{2})=1$ (whatever $i$) and
gets reward $r_{wait}$. With the action “search”, the robot goes from state
“high” to state “low” with probability $p(s_{1}|s_{2},a_{1})=1-\beta$ and gets
reward $r_{search}>r_{wait}$ with probability 1. In the state “low”, it stays
in the same state also with probability $\beta$ and gets reward $r_{search}$
in this case, but it gets the negative reward $r_{rescue}$ otherwise and goes
back to state “high”. Finally, with the action “recharge”, the robot goes to
state “high” with probability $p(s_{2}|s_{i},a_{3})=1$ and gets reward 0.

For the simulations, we consider the following parameters: $r_{rescue}=-1$,
$r_{wait}=0.4$, $r_{search}=0.9$, $\beta=0.8$, and the discount $\lambda=0.8$.
In this case, the maximum expected gain policy consists in performing the
action “search” for both states, i.e., the robot always searches, regardless
of the state. We also consider an alternative policy where the robot
implements the action “wait” in the state “low”.

With the alternative policy, the expected gain for state “low”
$\mathbb{E}[G_{t}|s_{1}]$ is reduced: 2 against 3.18. However,
$p(G_{t}>2-\epsilon|s_{1})$, $\epsilon>0$ being a small quantity, is greater
with the alternative policy: This can be observed on Figure 1 where we show
the empirical complement cumulative density function (CCDF) obtained by
simulating both policies. The maximum expected gain strategy is NOT optimal to
maximize $p(G_{t}>\alpha|s_{1})$ where $\alpha<2$.

Figure 1: Empirical CCDF, showing $p(G_{t}>x|s_{1})$, obtained with the
maximum expected gain policy and the alternative policy.

##  IV Proposed algorithm

We first introduce an algorithm to compute $v(s,\alpha)=p(G_{t}>\alpha|s)$
under a given policy (Algorithm 1). This algorithm is then slightly modified
to find a policy that maximizes $v(s,\alpha)$ (Algorithm 2).

###  IV-A Computing $v(s,\alpha)$ under a given policy

####  IV-A1 Notion of path to compute $v(s,\alpha)$

In this first subsection, we consider a given series of actions, e.g.,
established via the maximum expected gain policy. Figure 2 shows a trellis
with the transitions between the states at time $t$, $t+1$, and $t+2$, of an
arbitrary MDP.

Figure 2: Trellis representing the transitions between the states and the
rewards at time $t$, $t+1$, and $t+2$. The labels on the edges show the state
transition probabilities and short-term reward under several actions ($a_{1}$
and $a_{2}$).

We can define the notion of path in the trellis, corresponding to one possible
set of successive (environment) realization starting from a given state. The
$k$-th path starting from a state $s$ is defined by two values:

  * •

A path probability: $P_{s}(k)$.

  * •

A path gain: $G_{s}(k)$.

In the example of Figure 2, assume that the action $a_{1}$ is taken both at
time $t$ and $t+1$ (illustrated by the plain lines). The bold lines on the
figure show two paths: a path with $P_{s_{2}}(1)=p\cdot(1-p)$ and
$G_{s_{2}}(k)=-1+\lambda\cdot 0$, a path with $P_{s_{2}}(2)=p\cdot p$ and
$G_{s_{2}}(2)=-1+\lambda\cdot(-1)$.

Let us merge two paths (i.e., add their probability values) if their gain
difference is smaller than a some small value $\epsilon$. Consequently, in the
discounted infinite-horizon scope there is a finite number of paths with
distinct values, say $K$. We denote by

$\displaystyle P_{s}=[P_{s}(1),...,P_{s}(K)]\text{ and
}G_{s}=[G_{s}(1),...,G_{s}(K)]$ (12)

the vectors representing the probabilities and the gains, respectively, for
all paths starting from $s$. As a result, a state $s$ is characterized by the
set $\Omega_{s}=\\{P_{s},G_{s}\\}$. Then, $v(s,\alpha)$ can be computed as

$\displaystyle
v(s,\alpha)=\sum_{k=1}^{K}P_{s}(k)\cdot\mathbbm{1}\\{G_{s}(k)>\alpha\\},$ (13)

where $\mathbbm{1}\\{\cdot\\}$ denotes the indicator function.

####  IV-A2 Recursively computing $v(s,\alpha)$

Let us now consider one section of the trellis between a state $s_{1}$ and its
subsequent states (assuming an arbitrary action): $s_{1}^{\prime}$ with
probability $p$ and reward111For the sake of simplicity we write
$r(s_{1}^{\prime})$ for $r(s_{1},s^{\prime}_{1})$. $r(s_{1}^{\prime})$ and
$s_{2}^{\prime}$ with probability $1-p$ and reward $r(s_{2}^{\prime})$, as
illustrated on Figure 3. Assume also that
$\Omega_{s^{\prime}_{1}}=\\{P_{s^{\prime}_{1}},G_{s^{\prime}_{1}}\\}$ and
$\Omega_{s^{\prime}_{2}}=\\{P_{s^{\prime}_{2}},G_{s^{\prime}_{2}}\\}$ are
known and both vectors in $\Omega_{s^{\prime}_{1}}$ and
$\Omega_{s^{\prime}_{2}}$ are of size $K$.

Figure 3: One section of a trellis between a state $s_{1}$ and its subsequent
states.

The set $\Omega_{s_{1}}=\\{P_{s_{1}},G_{s_{1}}\\}$ can be computed from
$\Omega_{s^{\prime}_{1}}$ and $\Omega_{s^{\prime}_{2}}$ as follows:

$\displaystyle P_{s_{1}}=[p\cdot P_{s^{\prime}_{1}},(1-p)\cdot
P_{s^{\prime}_{2}})],$ (14) $\displaystyle
G_{s_{1}}=[r(s^{\prime}_{1})+\lambda\cdot
G_{s^{\prime}_{1}},r(s^{\prime}_{2})+\lambda\cdot G_{s_{2}^{\prime}}].$ (15)

The main drawback of the above equation222The equation also brings out the
main drawback of (12): as the size doubles when moving from one trellis
section to the previous one, $K$ is exponential in the trellis depth. is that
the size of the two resulting vectors is doubled compared to the ones of
$\Omega_{s^{\prime}_{1}}$ and $\Omega_{s^{\prime}_{2}}$.

We propose to adopt a binning strategy to maintain the size of the vectors
fixed. It is similar to the above rule where we merge two paths if their gain
difference is smaller than some value $\epsilon$. We introduce the vector
$G^{ref}_{s}$ whose components represent the center of the bins. The number of
bins $K$ (which is different from the number of paths $K$ in (12)), their
width, and their center is determined offline. For instance, the value of
$G^{ref}_{s_{1}}(K/2)$ can be chosen as $v(s_{1})$. Note that a unique vector
$G^{ref}_{s}$, $\forall s\in S$, could also be used to reduce the complexity.

Hence, after computing (14) and (15), if two components of $G_{s_{1}}$ fall in
the same bin of $G^{ref}_{s_{1}}$, e.g., the one centered at
$G^{ref}_{s_{1}}(k)$, the corresponding probability values are added and the
resulting value $P_{s_{1}}(k)$ is stored. To summarize, $G_{s_{1}}$ and
$G_{s_{1}}^{ref}$ are used to establish the binning rule, where $G_{s_{1}}$ is
itself computed from $G^{ref}_{s^{\prime}_{1}}$ and $G^{ref}_{s_{2}^{\prime}}$
as (instead of (15))

$\displaystyle G_{s_{1}}=[r(s^{\prime}_{1})+\lambda\cdot
G^{ref}_{s^{\prime}_{1}},r(s^{\prime}_{2})+\lambda\cdot
G^{ref}_{s_{2}^{\prime}}].$ (16)

As a result, we can use (14) and (16), as well as the binning trick, to
iteratively approximate $\Omega_{s}$. This is summarized within Algorithm 1.

Algorithm 1 Computing $\Omega_{s}$ and $v(s,\alpha)$, $\forall s\in S$, with a
given policy $\pi$.

Input: $G^{ref}_{s}$ with $K$ elements, $\forall s\in~{}S$, a policy $\pi$,
and $\alpha$.

1: Initialize $P_{s}$, $\forall s\in S$.

2: For all $s\in S$, choose $a$ according to $\pi$ and compute  $\displaystyle
G_{s}=[r(s_{j_{1}})+\lambda\cdot
G^{ref}_{s_{j_{1}}},...,r(s_{j_{k}})+\lambda\cdot G^{ref}_{s_{j_{k}}},...],$
(17) for all states $s_{j_{k}}$ with a non-zero transition probability
$p(s_{j_{k}}|s,a)$.

3: while a stopping criterion is not met do

4: for all states $s\in S$ do

5: Choose $a$ according to $\pi$ and compute  $\displaystyle
P_{s}=[p(s_{j_{1}}|s,a)\cdot P_{s_{j_{1}}},...,p(s_{j_{k}}|s,a)\cdot
P_{s_{j_{k}}},...],$ (18) for all states $s_{j_{k}}$ with a non-zero
transition probability $p(s_{j_{k}}|s,a)$.

6: Reduce the size of the vector $P_{s}$ to $K$ elements via the binning of
$G_{s}$ using $G_{s}^{ref}$: If two values of $G_{s}$ fall in the same bin,
merge the corresponding probabilities by adding them.

7: end for

8: end while

9: Compute $v(s,\alpha)$, $\forall s\in S$, using $\Omega_{s}$ via (13).

10: Return $\Omega_{s}$ and $v(s,\alpha)$, $\forall s\in S$.

Since the binning rule is constant, the same probabilities of $P_{s}$ are
merged at each iteration. Hence, steps 5 and 6 could be merged by directly
computing a vector $P^{bin}_{s}$ of size $K$ where the probability values are
added according to the binning rule.

Example of binning rule: On Figure 3, assume that
$G_{s_{1}}^{ref}=G_{s^{\prime}_{1}}^{ref}=G_{s^{\prime}_{2}}^{ref}=[1,\ 3]$,
$r(s^{\prime}_{1})=2$, $r(s^{\prime}_{2})=0$, and $\lambda=0.8$. Then, (16)
yields $G_{s_{1}}=[2.8,\ 4.4,\ 0.8,\ 2.4]$. Applying the binning rule yields
$P^{bin}_{s_{1}}=[P_{s_{1}}(3),\ P_{s_{1}}(1)+P_{s_{1}}(2)+P_{s_{1}}(4)].$

Note that $\Omega_{s}$ can also be used to compute the CCDF, as an alternative
to the empirical CCDF, without formerly running the algorithm.

###  IV-B Finding a policy optimized for $p(G_{t}<\alpha)$

In the previous subsection, the actions are chosen according to a given policy
$\pi$, such as the maximum expected gain policy. We now explain how to proceed
to learn actions maximizing $v(s,\alpha)$.

At steps 5-6 of Algorithm 1, $\Omega_{s}$ is computed for only one action $a$
determined by $\pi$. Alternatively, one could compute $\Omega_{s,a}$ for all
$a\in A_{s}$. With the obtained vectors at iteration $n$, one can compute

$\displaystyle Q_{n}(s,a,\alpha)=\sum_{k}P_{s,a}(k)\cdot
1\\{G^{ref}_{s}(k)>\alpha\\}.$ (19)

Then, the action can be chosen as

$\displaystyle a^{*}=\underset{a\in A_{s}}{\text{arg max }}Q_{n}(s,a,\alpha),$
(20)

and

$\displaystyle v_{n}(s,\alpha)=Q_{n}(s,a^{*},\alpha).$ (21)

As a result, Algorithm 1 is modified to result in Algorithm 2 to find a policy
optimized for $p(G_{t}>\alpha)$. Note that steps 5-6 of Algorithm 1 are merged
into a unique step 5 where the binning rule is directly applied to get the
vector $P_{s,a}^{bin}$, as described at the end of the previous subsection.

The ouput of Algorithm 2 (which can be seen as a training phase similarly to
the value iteration algorithm) is $\Omega_{s}$. This can then be used in the
inference phase to recover the optimal action to be performed at each state
(as $v(s)$ is used to get the action for the maximum expected gain policy, see
the end of Section II).

Algorithm 2 Finding a policy optimized for $p(G_{t}>\alpha)$.

Input: $\alpha$ and $G^{ref}_{s}$ with $K$ elements, $\forall s\in S$.

1: Initialize $P_{s}$, $\forall s\in S$, and set $n=0$.

2: For all $s\in S$ and for all $a\in A_{s}$ compute  $\displaystyle
G_{s,a}=[r(s_{j_{1}})+\lambda\cdot
G^{ref}_{s_{j_{1}}},...,r(s_{j_{k}})+\lambda\cdot G^{ref}_{s_{j_{k}}},...],$
(22) for all states $s_{j_{k}}$ with non-zero transition probabilities
$p(s_{j_{k}}|s,a)$ and establish the binning rules accordingly.

3: while a stopping criterion is not met do

4: for all states $s\in S$ do

5: For all $a\in A_{s}$ compute
$\displaystyle\begin{split}&P^{bin}_{s,a}=[p(s_{j_{1}}|s,a)P_{s_{j_{1}}}(k_{1})+p(s_{j_{2}}|s,a)P_{s_{j_{2}}}(k_{2})+\\\
&...,...].\end{split}$ (23) for all states $s_{j_{k}}$ with non-zero
transition probabilities $p(s_{j_{k}}|s,a)$ and where the probabilities are
added according to the binning rule of step 2.

6: For all $a\in A_{s}$ compute  $\displaystyle
Q^{n}(s,a,\alpha)=\sum_{k}P^{bin}_{s,a}(k)\cdot 1\\{G^{ref}_{s}(k)>\alpha\\},$
(24) and set $a^{*}=\underset{a\in A_{s}}{\text{arg max }}Q(s,a,\alpha)$.

7: Set $P_{s}=P^{bin}_{s,a^{*}}$.

8: end for

9: Increment $n$.

10: end while

11: Return Return $\Omega_{s}$, $\forall s\in S$.

###  IV-C Simulation results on the toy example

Via Figure 1 for the recycling robot example, we observe that the optimal
policy for $p(G_{t}>\alpha|s_{1})$ changes if $\alpha$ is greater or smaller
than 2. As an example, we run Algorithm 2 with $\alpha=1.8$ and $\alpha=2.2$.
In the first case, the CCDF of the policy obtained is the one shown with the
blue curve on Figure 1. In the second case, we get the red curve. This yields
the expected results, meaning that Algorithm 2 manages to find the optimal
strategy.

##  V Towards a deep learning extension

###  V-A Temporal-difference learning

A famous alternative to the value iteration algorithm is temporal-difference
learning, see [5, Chap. 6], where the model is updated based on the difference
between two estimates of a state. As an example, the $Q$ learning algorithm
relies on this paradigm. Temporal-difference learning also contains the main
idea behind deep reinforcement learning algorithms, such as deep $Q$ learning
[2], as it enables to produce an error signal to train the neural networks.

Instead of looping over all the states to estimate the $v^{\pi^{*}}(s)$, the
algorithm walks from state to state. In the standard implementation using the
Monte Carlo approach, the system is in a state $s$ at time $t$ and an action
$a$ is chosen by the algorithm. One then gets a reward $r(s,s^{\prime})$ and
the new state $s^{\prime}$ of the system at time $t+1$. The estimate
$\hat{Q}(s,a)$ of $Q^{\pi^{*}}(s,a)$ is updated as

$\displaystyle\hat{Q}(s,a)=\hat{Q}(s,a)+\gamma\Delta_{t},$ (25)

instead of (8) with the value iteration algorithm. The quantity $\gamma$ is
the learning rate and

$\displaystyle\Delta_{t}=\hat{Q}(s,a)-\hat{Q}^{\prime}(s,a)$ (26)

is the error signal, where:

  * •

$\hat{Q}^{\prime}(s,a)=r(s,s^{\prime})+\lambda\max_{a^{\prime}}\hat{Q}(s^{\prime},a^{\prime})$
is a first estimate of $Q^{\pi^{*}}(s,a)$ based on the observed reward and
subsequent state obtained when taking the action $a$ in state $s$.

  * •

$\hat{Q}(s,a)$ is a second estimate of $Q^{\pi^{*}}(s,a)$ obtained via the
currently stored value.

In our framework, a similar error signal can be generated with the vector of
probabilities. We use the temporal-difference idea but without the Monte Carlo
approach as we assume knowledge of the model. For a state $s$ and given a
binning rule, $P_{s}$ can be updated as

$\displaystyle\begin{split}&P_{s}=P_{s}+\gamma\Delta_{t},\end{split}$ (27)

where $\Delta_{t}=f(P^{bin}_{s},P_{s})$ is the error signal with:

  * •

$P^{bin}_{s}$ is obtained via (23) and (24),

  * •

$P_{s}$ is the currently stored value.

The function $f$ could be a distance measure between two distributions such as
the KL divergence or the squared norm of the difference of the two vectors. Of
course, (27) needs to be normalized at each time step.

The drawback here, compared to standard temporal-difference learning relying
on Monte Carlo method (and thus “true” reinforcement learning), is that we
still need a model of the environment dynamic (i.e., we need the transition
probablities to compute (23)). For an algorithm compliant with the
reinforcement learning paradigm, see Algorithm 1 in [1] and Appendix VII.
Alternatively, a neural network could be trained in a Monte Carlo phase to
infer the transition probabilities and then used within the algorithm
described above.

###  V-B With neural networks

The main idea behind deep $Q$ learning consists in using a neural network to
compute $v(s)$ or $Q(s,a)$, $\forall s\in S$, and using the error signal
$\Delta_{t}$ (or a sum of several error signals) to train the neural network.

We can propose a similar approach where a neural network computes $P_{s}$ and
the error signal $\Delta_{t}=f(P^{bin}_{s},P_{s})$ is used for the training.

##  VI Conclusions

In this paper, we introduced an algorithm to find a policy that maximizes the
probability $p(G_{t}>\alpha)$ that the gain $G_{t}$ is higher than a
predefined value $\alpha$. This is an alternative to algorithms searching for
a policy that maximizes the expected gain. Optimizing with respect to this new
criterion $p(G_{t}>\alpha)$ induces computing path gains and probabilities,
where a path corresponds to a given series of state transitions. Computing
such path metrics is intractable as the number of paths is exponential in the
depth of the trellis. As a result, we introduced a recursive calculation
method and a binning rule to merge paths having similar gains. While this new
algorithm is presented as an extension of the value iteration algorithm, the
main principles are not restricted to this paradigm. We show in the last
section how these principles can be applied to modify the $Q$ learning and
deep $Q$ learning algorithms to optimize the alternative objective. We use
temporal-difference learning but without the Monte Carlo approach.

##  VII Appendix

We discuss why the Monte Carlo approach seems to be efficient, as reported in
[1].

Consider the example of Figure 3, with rewards
$r(s^{\prime}_{1})=r(s^{\prime}_{2})=0$. Let $P_{s^{\prime}_{1}}$ and
$P_{s^{\prime}_{2}}$ be the vector of probabilities at $s^{\prime}_{1}$ and
$s^{\prime}_{2}$ respectively. Let $P^{*}_{s_{1}}=p\cdot
P_{s^{\prime}_{1}}+(1-p)\cdot P_{s^{\prime}_{2}}$ be the optimal probability
vector to learn at $s_{1}$ via Monte Carlo samples and $P^{model}_{s_{1}}$ an
estimation by a model.

In the Monte Carlo approach, the error signal is computed based on a
realization $(s,a,s^{\prime}_{i})$. If $s^{\prime}=s^{\prime}_{1}$, the error
signal is $f(P^{model}_{s_{1}},P_{s^{\prime}_{1}})$, where
$P_{s^{\prime}_{1}}$ acts as the label. If $s^{\prime}=s^{\prime}_{2}$, the
error signal is $f(P^{model}_{s_{1}},P_{s^{\prime}_{2}})$, where
$P_{s^{\prime}_{2}}$ acts as the label. In the Monte Carlo training process,
the learning algorithm gets approximately $p$ times the error
$f(P^{model}_{s_{1}},P_{s^{\prime}_{1}})$ and $(1-p)$ times the error
$f(P^{model}_{s_{1}},P_{s^{\prime}_{2}})$. Does it make $P^{model}_{s_{1}}$
converge to $P^{*}_{s_{1}}$?

We performed simple simulations, where $P_{s^{\prime}_{1}}$ and
$P_{s^{\prime}_{2}}$ are chosen as discrete Gaussians with distinct means and
where we train a neural network as follows. We use the gradient of
$f(P^{model}_{s_{1}},P_{s^{\prime}_{i}})$ to update the model, based on a
dataset containing $p$ labels $P_{s^{\prime}_{1}}$ and $(1-p)$ label
$P_{s^{\prime}_{2}}$ (and where the input of the neural network is a
constant). The output of the model $P^{model}_{s_{1}}$ converges towards
$P^{*}_{s_{1}}$, which explains why the Monte Carlo approach is also
efficient. There may be however some pathological cases.

## References

  * [1] M. Bellemare, W. Dabney, and R. Munos, “A distributional perspective on reinforcement learning,” International conference on machine learning, pp. 449-458, July 2017.
  * [2] Mnih et al., “Human-level control through deep reinforcement learning,” Nature, vol. 518, no. 7540, pp. 529–533, 2015.
  * [3] M. L. Puterman, “Markov Decision processes: Discrete Stochastic Dynamic Programming,” John Wiley $\&$ Sons, 1994.
  * [4] J.-C. Sibel, N. Gresset, and V. Corlay, “An Application-oriented scheduler,” 2023 IEEE Wireless Communications and Networking Conference (WCNC), March 2023. Online: https://arxiv.org/abs/2302.09926
  * [5] R. S. Sutton and A. G. Barto, “Reinforcement Learning, An Introduction”, 2nd ed., The MIT Press, 2020.


# A Quantum-Inspired Classical Solver for Boolean $k$-Satisfiability Problems

S. Andrew Lanham, Brian R. La Cour  Applied Research Laboratories
The University of Texas at Austin
Austin, TX, 78713-8029 USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we detail a classical algorithmic approach to the
$k$-satisfiability ($k$-SAT) problem that is inspired by the quantum amplitude
amplification algorithm. This work falls under the emerging field of quantum-
inspired classical algorithms. To propose our modification, we adopt an
existing problem model for $k$-SAT known as Universal SAT (UniSAT), which
casts the Boolean satisfiability problem as a non-convex global optimization
over a real-valued space. The quantum-inspired modification to UniSAT is to
apply a conditioning operation to the objective function that has the effect
of “amplifying” the function value at points corresponding to optimal
solutions. We describe the algorithm for achieving this amplification, termed
“AmplifySAT,” which follows a familiar two-step process of applying an oracle-
like operation followed by a reflection about the average. We then discuss
opportunities for meaningfully leveraging this processing in a classical
digital or analog computing setting, attempting to identify the strengths and
limitations of AmplifySAT in the context of existing non-convex optimization
strategies like simulated annealing and gradient descent.

©2021 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

##  I Introduction

The development of new quantum algorithms has given rise to rich tradeoffs
between classical and quantum solutions for different problems. Quantum
computing’s promise lies in the fact that there are many problems for which
quantum algorithms produce a significant speedup over the best known classical
algorithms. Yet at the same time, careful scrutiny of the differences between
quantum and classical algorithms has led to new classical analogues of quantum
algorithms, sometimes promising the same computational advantages. The
resulting area of research, termed quantum-inspired algorithms, means that our
understanding of the delineation between quantum and classical algorithms is
continually evolving. Quantum-inspired design approaches have produced novel
algorithms in a variety of application areas, perhaps most concretely in
recommendation systems [1], but also in reinforcement learning [2], constraint
satisfaction [3], and more general areas of machine learning and computational
linear algebra [4, 5, 6, 7].

Recent quantum-inspired algorithms leverage similarities between quantum state
preparation assumptions and input assumptions for classical algorithms [8, 4].
In these works, quantum state preparation assumptions are shown to be
equivalent to a set of classical assumptions involving sampling access to a
vector or matrix. Sampling access refers to a situation in which a vector or
matrix can be normalized and the entries treated as a discrete distribution to
be sampled from. The similarity between sampling access and state preparation
allows for a “dequantization” of a variety of quantum algorithms by adopting
the assumptions in a classical setting and proposing analogous classical
processing. A thorough characterization of this framework with applications to
quantum and classical machine learning algorithms is contained in [4].

Other families of quantum algorithms are much less likely to admit a quantum-
inspired classical analog because their speedup relies on uniquely quantum
resources like entanglement. In this work, we focus on an algorithm inspired
by amplitude amplification (AA). At a high level, AA uses principles of
quantum coherence to evolve a quantum state. Coherence is not a uniquely
quantum phenomenon, and, as we will explore in this work, there is potential
to use it to provide speedups for classical solvers. Additionally, the focus
on AA places this work in a different area than recent dequantized algorithms
based on quantum state preparation assumptions.

AA algorithms find application in areas including unstructured search [9],
function optimization [10], and some types of numerical estimation [11, 12].
More generally, AA can be applied to the output state of almost any quantum
algorithm to boost the probability of success for that algorithm [13]. Each
algorithm in the AA family works by increasing the amplitude of a desirable
element in a quantum superposition by $\mathcal{O}(1/\sqrt{N})$ for each loop
of the algorithm. It follows that with $\mathcal{O}(\sqrt{N})$ repetitions of
this algorithm, the desired element is measured with probability
$\mathcal{O}(1)$. In this sense, AA is a way of boosting the probability of
measuring a given element in a set (typically large) of possibilities, and as
such, is intimately related to sampling access to a vector.

In a classical setting, therefore, it is not immediately obvious how to
introduce an algorithm that is designed for sampling access. It is rare that
sampling access to classical data occurs naturally, and thus it must often be
considered as a pre-processing assumption, as occurs in much quantum-inspired
work [8, 4]. In contrast, classical data can often be assumed to admit
efficient query access, which simply refers to the fact that individual
entries of a vector or matrix can be efficiently retrieved. Since this work
details a classical analog to a quantum AA algorithm, we are forced to address
the presence or lack of sampling access. In this work we make no assumption of
sampling access because it would be far too costly and unnatural for the
problem setting. This leads to an exploration of the interplay between query
access and sampling access; in particular, we explore whether algorithms
designed assuming sampling access can be meaningfully used in situations in
which only query access is available. A quantum-inspired AA algorithm is a
prime candidate for exploring this question.

Boolean satisfiability is an important class of problem that can be treated as
an unstructured or semi-structured search, and is therefore amenable to AA.
Boolean satisfiability concerns questions of determining whether there exists
a binary string $x\in\\{0,1\\}^{n}$ for which a given Boolean function is
satisfied, that is, evaluates to $1$. Considerable effort has been focused on
developing better SAT solvers. For example, SAT solvers have been designed
based on a variety of heuristic search algorithms including backtracking [14,
15], conflict-driven clause learning (an evolution of backtracking) [16],
local search [17, 18, 19], and simulated annealing [20]. A classical
randomized algorithm for $3$-SAT problems achieving average complexity within
a polynomial factor of $\mathcal{O}(1.334^{n})$ was proposed by Schöning [21].
These algorithms and similar ones have also led to a number of hardware-based
SAT solvers achieving significant speedups [22, 23, 24]. With quantum
computers, Grover’s algorithm stands out as the most relevant algorithm for
solving SAT problems [9]. Its generalization, AA, can be used to achieve
quadratic speedups over the best known classical algorithms for search [13].

A separate algorithmic paradigm for solving SAT problems casts discrete
Boolean SAT problems as global minimizations over a bounded real-valued space
[25, 26, 27, 28]. With the conversion to a continuous domain, the resulting
minimization task becomes amenable to a variety of solution strategies from
optimization, such as gradient descent and its generalizations. Nevertheless,
the objective function that results from this conversion is likely non-convex
and so solution strategies remain largely heuristic in design. A SAT solution
strategy broadly based on the problem formulation of [29], known as UniSAT,
was proposed in [30] and an analog solver was presented in [22]. The work of
[30] significantly extends the original formulation by proposing Lagrange
multiplier methods as a strategy for avoiding local minima. A continuous-time
dynamical system formulation is presented with efficient convergence to the
global minimum, but at exponential energy cost. The work of [22] follows this
up with an analog circuit demonstration implementing this dynamical system. In
this work, the discussion of quantum-inspired algorithms with possible
applications to analog solvers loosely relates to the field of quantum
emulation, which proposes classical analog circuits for executing quantum
algorithms [31, 32].

This work continues development under the UniSAT problem model by proposing a
modification of the objective functions originally presented in [29] and [30].
Modifications to the objective function are common in optimization, for
example, when introducing barrier functions to convert a constrained
optimization to an unconstrained one. Our modification, termed AmplifySAT, is
inspired by the quantum AA algorithm, and has the qualitative effect of
“amplifying”, or “deepening” the global optimum of an objective function
encoding a SAT problem. Transformations of the objective function that target
the extrema are not common, but have been explored in some works, e.g., [33,
34]. The iterative nature of AA means that the degree of amplification of the
extrema can be coarsely controlled. We continue in Section II with a review of
the UniSAT problem formulation. In Section III we propose using AmplifySAT to
modify the objective function. We conclude in Section IV with an exploration
of the numerous possibilities to follow-up this investigation. In particular,
we suggest reconsidering the performance of various standard non-convex
optimization strategies with objective functions conditioned by AmplifySAT. In
this sense, this work does not propose a brand new optimization algorithm, but
simply views AmplifySAT as a conditioning stage that may improve the
performance of existing optimization approaches.

##  II The UniSAT Problem Model

We begin by reviewing the $k$-SAT problem and its UniSAT formulation. The goal
when solving such a problem is to find the values of $n$ variables
$\mathbf{x}\in\\{0,1\\}^{n}$ that satisfy a propositional formula $F$. The
formula is assumed to be in the conjunctive normal form, that is, a
conjunction of clauses $C_{m},m=1,\ldots,M$ with

$F(\mathbf{x})=\bigwedge_{m=1}^{M}C_{m}$ (1)

We have $F(\mathbf{x}):\mathbb{F}_{2}^{n}\mapsto\\{0,1\\}$, and the function
returns $1=\textrm{TRUE}$ if the function is satisfied for input $\mathbf{x}$,
and otherwise it returns $0=\textrm{FALSE}$. A clause is a disjunction between
$k$ variables. Complements of variables, e.g. $\bar{x}_{i}$ for variable
$x_{i}$ are valid variables as well. As an example, a $k$-SAT clause could
take the form

$C_{1}=(x_{1}\lor x_{4}\lor\bar{x}_{7})\,.$ (2)

The UniSAT problem model expresses a $k$-satisfiability problem as a function
over a bounded continuous space so that it can be solved using different
classes of solvers [29, 30]. In the $k$-UniSAT relaxation, each clause $C_{m}$
is represented by a continuous function
$K_{m}(\mathbf{s}):[-1,1]^{n}\mapsto[0,1]$ given by

$K_{m}(\mathbf{s})=2^{-k}\prod_{i=1}^{n}(1-c_{m,i}s_{i})$ (3)

The objective function is defined over the $n$-dimensional hypercube, i.e.
$[-1,1]^{n}$. The vertices of this hypercube correspond to the binary inputs
to the Boolean satisfiability problem of Equation (1). This means that
$\mathbf{s}$ obeys the equivalence that $s_{i}=-1\mapsto\textrm{FALSE}$ and
$s_{i}=+1\mapsto\textrm{TRUE}$ for each entry $s_{i}$ of $\mathbf{s}$. We
express the set of vertices for an $n$-variable problem as

$\mathcal{V}=\\{\mathbf{x}:\mathbf{x}\in\mathbb{R}^{n},x_{i}\in\\{-1,1\\}\;\forall
i\in[n]\\}\,.$ (4)

The clause variables $c_{m,i}$ take value $1$ when variable $x_{i}$ appears in
clause $m$, $-1$ if $\bar{x}_{i}$ appears, and $0$ if neither $x_{i}$ nor
$\bar{x}_{i}$ appear in clause $m$. From the clause functions
$K_{m}(\mathbf{s})$, [29] defines an objective function

$E(\mathbf{s})=\sum_{m=1}^{M}K_{m}(\mathbf{s})\,,$ (5)

with the property that if $E(\mathbf{s}^{*})=0$ for
$\mathbf{s}^{*}\in\mathcal{V}$ then $\mathbf{s}^{*}$ corresponds to a
satisfying solution. The satisfying solutions, therefore, are encoded as
vectors $\mathbf{s}\in\mathcal{S}^{*}\subseteq\mathcal{V}$ achieving the
global minimum. Here we have defined $\mathcal{S}^{*}$ as the set of vectors
mapping to satisfying solutions to a given satisfiability problem $F$.
Clearly, the objective function $E(\mathbf{s})$ may also be evaluated at
points in the hypercube that are not vertices to assist in finding the global
minimum.

##  III A Modification to the Objective Function Inspired by Amplitude
Amplification

The role of AA in satisfiability is to boost the probability amplitudes
associated with the satisfying solutions to a problem of interest. When
considering a quantum-inspired version of this algorithm, we must contend with
the numerous differences between classical and quantum processing.
Fundamentally the difference can be understood in terms of how information
encoded in a quantum or classical state is efficiently accessed. The types of
access for each setting are best understood as “sampling access” and “query
access,” and their relevance has been highlighted in other quantum-inspired
classical algorithms [8, 4]. We review these types of access now. The quantum
paradigm typically admits efficient sampling access to the set of feasible
solutions when they are encoded in a quantum state. This is because in a
quantum state, the formalism of probability amplitudes and quantum
measurements is integral to how data is accessed. Sampling access is
efficiently realized through the process of quantum measurement. With
measurement, each possible outcome is observed with frequency proportional to
the modulus-squared of the probability amplitude of that outcome. In a
classical setting there is no efficiently realizable analogy to this notion,
in general. In contrast, only the output values corresponding to individual
inputs can be “queried” efficiently, meaning that an individual input must be
specified and information (e.g. “satisfying” or “not satisfying”) for only
that input is returned.

The AmplifySAT algorithm presented in this section, being inspired by AA, is
naturally optimized for sampling access. Since there is no well-defined notion
of sampling from a classical objective function, we will have to address
whether an algorithm designed for sampling access is even useful in a setting
where only query access is efficiently available. Despite this apparent
incompatibility, we will see that AA has observable effects on query accesses
as well, and that it has the potential to improve the performance of
algorithms relying only on query access. In particular, the amplification step
introduced in AmplifySAT affects properties of the objective function relevant
to the convergence of non-convex solvers that only access the function
locally. The algorithm works as follows. First, we construct a new objective
function $V(\mathbf{s})$ based on the original objective function
$E(\mathbf{s})$, but for which the satisfying arguments
$\mathbf{s}\in\mathcal{S}^{*}$ are maxima instead of minima. AmplifySAT then
iteratively transforms $V(\mathbf{s})$ into a conditioned function
$V_{\ell}(\mathbf{s})$ for the $\ell^{\textrm{th}}$ iterate, whose extremal
points are “amplified” using a two-step process similar to multiplication by
an oracle operator followed by reflection about the average.

###  III-A Modified Objective Function

The objective function $E(\mathbf{s})$ of Equation (5) has a variety of useful
properties; namely, it’s constrained to the interval $[0,M]$ so it’s non-
negative. Furthermore, the set of binary inputs to the SAT problem can be
mapped to the vertices $\mathcal{V}$ of the $n$-dimensional hypercube
$[-1,1]^{n}$. But in order to benefit from amplification effects, we must
consider an objective function for which the optimal points are the global
maxima. This is easily obtained as a modification of $E(\mathbf{s})$ by
substituting each clause function $K_{m}(\mathbf{s})$ by
$1-K_{m}(\mathbf{s})$, giving

$V(\mathbf{s})=\sum_{m=1}^{M}(1-K_{m}(\mathbf{s}))\,.$ (6)

Since each $K_{m}(\mathbf{s})\in[0,1]$, each term $(1-K_{m}(\mathbf{s}))$ in
the summation is also in the same interval. The minima of $E(\mathbf{s})$ now
appear as maxima in $V(\mathbf{s})$. When all clauses are fully satisfied at a
vertex, we have $V(\mathbf{s})=M$, and the global maximum is attained here,
i.e. $\mathbf{s}\in\mathcal{S}^{*}$.

The value of $V(\mathbf{s})$ at the points corresponding to binary solutions
$\mathbf{s}\in\mathcal{V}$ indicates the cost for that binary input. A list of
these costs could theoretically be expressed in a vector
$\mathbf{c}\in\mathbb{R}^{2^{n}}$. The resulting vector $\mathbf{c}$ is
loosely similar to a quantum state vector encoding the problem, in the sense
that we could normalize and sample from it, and each entry could be mapped to
the set of classical solutions. In this case, the modes of the resulting
distribution would correspond to the optimal solutions, since the objective
function attains maximum values at those points. It was shown in [13] that AA
can be performed directly on the output state of an arbitrary quantum
algorithm in order to boost the probabilities corresponding to valid
solutions. The initial output state could be a superposition of all possible
outcomes, such as with Grover’s algorithm. The same concept turns out to apply
to un-normalized vectors as well. In summary, if we had direct sample access
to a vector like $\mathbf{c}$, constructed by sampling the value of the
objective function at $2^{n}$ points, we could apply Grover-inspired oracle
and diffusion operators to that vector to boost the amplitudes of entries
corresponding to solutions only. In the classical setting we don’t have access
to such a vector, except through single queries to the objective function
$V(\mathbf{s})$ for individual $\mathbf{s}\in\mathcal{V}$. But we can still
define a set of transformations to the objective function such that each query
for $\mathbf{s}\in\mathcal{V}$ will return a conditioned cost
$V_{\ell}(\mathbf{s})$ according to an AA algorithm.

Algorithm 1 AmplifySAT($\mathbf{s},\ell$)

0: $\mathbf{s}\in[-1,1]^{n}$

0: $y\in\mathbb{R}$

1: $y\leftarrow f_{0}(\mathbf{s})$

2: if ($\ell>0$) then

3: for $i=0$ to $\ell$ do

4: $z\leftarrow 2\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})y]-T(\mathbf{s})y$

5: $y\leftarrow z$

6: end for

7: end if

8: return $y$

###  III-B Amplitude-Amplification-Inspired Conditioning

The steps of AmplifySAT closely resemble its quantum analogs. In particular,
we first propose multiplying $V(\mathbf{s})$ by a function $T(\mathbf{s})$
whose action on the vertices of the hypercube is identical to the action of a
quantum oracle operator $U_{\omega}$ on a quantum state
$|\psi\rangle\in\mathbb{C}^{2^{n}}$. We then follow this by a reflection about
the average value of the function. A valid oracle function, $T(\mathbf{s})$,
must take the following values at the vertices $\mathbf{s}\in\mathcal{V}$:

$T(\mathbf{s})=\begin{cases}-1\quad&\mathbf{s}\in\mathcal{S}^{*}\\\
1\quad&\mathbf{s}\in\mathcal{V}\setminus\mathcal{S}^{*}\,.\end{cases}$ (7)

The action of $T(\mathbf{s})$, then, when used to form
$T(\mathbf{s})V(\mathbf{s})$ is to flip the sign of the objective function at
points where $\mathbf{s}$ corresponds to a satisfying classical input, i.e.
$\mathbf{s}\in\mathcal{S}^{*}$, the set of satisfying inputs, and to do
nothing to the function value at an unsatisfying argument.

There are many ways to construct a function like $T(\mathbf{s})$. One such
construction follows simply from the clause functions $K_{m}(\mathbf{s})$. In
particular, we recognize that $K_{m}(\mathbf{s})=0$ for a fully satisfied
clause and $K_{m}(\mathbf{s})=1$ for a fully unsatisfied clause. Then
trivially we have the following relation for the product of clause functions:

$\prod_{m=1}^{M}(1-K_{m}(\mathbf{s}))=\begin{cases}1&\mathbf{s}\in\mathcal{S}^{*}\\\
0&\mathbf{s}\in\mathcal{V}\setminus\mathcal{S}^{*}\end{cases}\,.$ (8)

The product of clause functions enforces the conjunction between clauses,
since one or more unsatisfying clauses will render a $K_{m}(\mathbf{s})$ equal
to $1$, and thus the entire product of Equation (8) equal to $0$, or
unsatisfied. From now on, we make the substitution
$\tilde{K}_{m}(\mathbf{s})=1-K_{m}(\mathbf{s})$. Equation (8) can now be
shifted and scaled to form the phase-flip-style oracle function
$T(\mathbf{s})$ as

$T(\mathbf{s})=-2(\prod_{m=1}^{M}\tilde{K}_{m}(\mathbf{s})-\frac{1}{2})$ (9)

and $T(\mathbf{s})$ satisfies the description of Equation (7). Of course, the
oracle function is nonlinear in the variables, and may have many local extrema
within the interior of the hypercube.

Other formulations of oracle functions are possible, so long as the result
behaves like Equation (7) at the vertices. Another satisfying choice is

$T_{2}(\mathbf{s})=-2(\min_{m\in[M]}\tilde{K}_{m}(\mathbf{s})-\frac{1}{2})\,.$
(10)

Each unsatisfying solution $\mathbf{s}\in\mathcal{V}\setminus\mathcal{S}^{*}$
will have a fully unsatisfied clause, rendering $\tilde{K}_{m}(\mathbf{s})=0$
for some $m$. In contrast, at the satisfying solution we have
$\min_{m\in[M]}\tilde{K}_{m}(\mathbf{s})=1$ and $T_{1}(\mathbf{s}^{*})=-1$.
This meets the criteria of Equation (7). Another example that might be used
for its gradient characteristics is

$T_{3}(\mathbf{s})=-2(\frac{V(\mathbf{s})}{M}\min_{m\in[M]}\tilde{K}_{m}(\mathbf{s})-\frac{1}{2})\,.$
(11)

We expect that each choice of oracle function modifies the optimization
landscape differently in the full algorithm, affecting properties related to
the gradient, smoothness, number of local minima, and thus the overall
efficacy of solvers.

After pointwise multiplication by $T(\mathbf{s})$, AmplifySAT continues
mirroring AA algorithms by reflecting the function about its average value.
But, before we detail this second stage, we also remark on the initial
function by which to pointwise multiply $T(\mathbf{s})$. Throughout this
section, we have suggested starting with objective function $V(\mathbf{s})$.
But as previously mentioned, the effects of amplitude amplification
essentially work on the output of any quantum algorithm $\mathcal{A}$ [13].
That is, an arbitrary quantum algorithm $\mathcal{A}$ producing a
superposition of classical states
$\mathcal{A}|0\rangle=\sum_{i}a_{i}|i\rangle$ can be subjected to AA to boost
the probability of measuring a “good” solution. In Grover’s algorithm,
$\mathcal{A}=\mathbf{H}_{2^{n}}$, with the interpretation that the Hadamard
matrix is an algorithm producing a uniform distribution over all possible
solutions.

In the analogous quantum-inspired setting, therefore, we might also consider
different choices for the initial objective function, which we now refer to as
$f_{0}(\mathbf{s})$, that reflect a state of prior knowledge. For example, the
choice $f_{0}(\mathbf{s})=1$ represents an initial function that reflects no
structure to the search problem and resulting optimization landscape, and is
thus analogous to $\mathcal{A}=\mathbf{H}_{2^{n}}$. In this case, repeated
application of the oracle function $T(\mathbf{s})$ and reflection about the
average would do all the work of encoding information about the optimal
solutions. We may alternatively proceed with our choice
$f_{0}(\mathbf{s})=V(\mathbf{s})$ as was suggested by the initial example,
since it is an objective function that encodes useful initial information
about the SAT problem; in particular, the cost of this function at each
$\mathbf{s}\subseteq\mathcal{V}$ reports the number of clauses satisfied by
that input. This introduces an intermediate measure of value of each solution,
where points that have a higher value $V(\mathbf{s})$ are “closer” to
optimal.” However, as we will see, it is also necessary in AmplifySAT to pre-
compute the average value of the initial function $f_{0}(\mathbf{s})$, which
can limit the available options. Fortunately, the choices
$f_{0}(\mathbf{s})=1$ and $f_{0}(\mathbf{s})=V(\mathbf{s})$ both admit pre-
computed averages.

###  III-C Averages

In the previous discussion we see a melding of the ideas of a classical
objective function and the $2^{n}$-ary outcomes that are encoded into a
quantum state vector. The action of the oracle function $T(\mathbf{s})$ is
identical to the action of the quantum oracle operator $U_{\omega}$ for a
theoretical vector of costs. It turns out we can also apply a reflection about
the average in our inspired algorithm as well. In fact, the average values are
directly derivable for each iteration of the algorithm offline and in closed
form, as long as the average values of $f_{0}(\mathbf{s})$ and $T(\mathbf{s})$
at the vertices are also known. For the basic construction of the algorithm,
we continue to be primarily concerned with the behavior of the functions at
the vertices only, because these are points corresponding to binary solutions.
For this reason it is sufficient to define the average value of a function as

$\mathbb{E}_{\mathbf{s}}[f(\mathbf{s})]=\frac{1}{2^{n}}\sum_{\mathbf{s}\in\mathcal{V}}f(\mathbf{s})\,.$
(12)

as opposed to a cumbersome integral over the entire interior $[-1,1]^{n}$. We
use the expectation operator for notational convenience; there are no random
variables arising in Equation (12).

We may now define $s(\cdot):\mathbb{R}\mapsto\mathbb{R}$ as the linear
functional performing reflection about the average. That is, in general,

$s\circ
f(\mathbf{s})=2\mathbb{E}_{\mathbf{s}}[f(\mathbf{s})]-f(\mathbf{s})\,.$ (13)

One full iteration of AmplifySAT would thus produce the function

$s\circ(T(\mathbf{s})f_{0}(\mathbf{s}))=2\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{0}(\mathbf{s})]-T(\mathbf{s})f_{0}(\mathbf{s})$
(14)

In a classical setting, we cannot apply such a reflection unless we know the
value of the constant
$\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{0}(\mathbf{s})]$ in advance, or we’d
have to sample the function $T(\mathbf{s})f(\mathbf{s})$ to estimate the
average value. Fortunately, as we now show, this average is available in
closed form for the initial functions $f_{0}(\mathbf{s})$ of greatest
interest, as long as the number of variables $n$, number of satisfying
solutions $L$, number of clauses $M$, and averages
$\mathbb{E}_{\mathbf{s}}[f_{0}(\mathbf{s})]$ and
$\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})]$ are known. Such prior assumptions are
typical of AA algorithms and of satisfiability problems. Also, when the number
of solutions $L$ is not known, there are still methods for applying AA
algorithms effectively [13].

We can pre-compute $\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{0}(\mathbf{s})]$
for these problems by making the following observations. First, we note that
the action of pointwise multiplying by $T(\mathbf{s})$ will flip the sign of
the amplitude at the satisfying solutions. The value of $f_{0}(\mathbf{s})$ at
these satisfying solutions is originally $f_{0}(\mathbf{s}^{*})$. Thus the
average value of the function with the flipped sign should change to

$\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{0}(\mathbf{s})]=\mathbb{E}_{\mathbf{s}}[f_{0}(\mathbf{s})]-2\frac{f_{0}(\mathbf{s}^{*})L}{2^{n}}\,.$
(15)

We subtract one term of $f_{0}(\mathbf{s}^{*})L/2^{n}$ to remove the
contributions of the positive-sign optimal terms in
$\mathbb{E}_{\mathbf{s}}[f_{0}(\mathbf{s})]$, and then we subtract the second
term to reflect that the sign at those points is now negative (i.e. amplitudes
of $-1$). This gives us a closed form expression of
$\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{0}(\mathbf{s})]$ for many choices of
$f_{0}(\mathbf{s})$. For instance, with our choice of $f_{0}(\mathbf{s})=1$ we
have $\mathbb{E}_{\mathbf{s}}[f_{0}(\mathbf{s})]=1$ and
$f_{0}(\mathbf{s}^{*})=1$.

For the case $f_{0}(\mathbf{s})=V(\mathbf{s})$ we can easily deduce that
$f_{0}(\mathbf{s}^{*})=M$ and also calculate that

$\displaystyle\mathbb{E}_{\mathbf{s}}[V(\mathbf{s})]$
$\displaystyle=\mathbb{E}_{\mathbf{s}}\biggl{[}\sum_{m=1}^{M}\bigl{(}1-2^{-k}\prod_{i=1}^{n}(1-c_{m,i}s_{i})\bigr{)}\biggr{]}$
(16a)
$\displaystyle=\sum_{m=1}^{M}\mathbb{E}_{\mathbf{s}}\biggl{[}\bigl{(}1-2^{-k}\prod_{i=1}^{n}(1-c_{m,i}s_{i})\bigr{)}\biggr{]}$
(16b)
$\displaystyle=M-\sum_{m=1}^{M}2^{-k}\mathbb{E}_{\mathbf{s}}[\prod_{i=1}^{n}(1-c_{m,i}s_{i})]\,.$
(16c)

Recall that the product
$K_{m}(\mathbf{s})=2^{-k}\prod_{i=1}^{n}(1-c_{m,i}s_{i})$ equals zero when the
clause is satisfied and $1$ when unsatisfied. Any $k$ variables or their
negations featured in a clause will be unsatisfied for exactly $2^{n-k}$
cases, since for any disjunction of $k$ variables, there is only one
combination of those variables that makes that disjunction equal to $0$. That
leaves $2^{n}-2^{n-k}$ cases for which the clause is satisfied. In
expectation, therefore, we have for any clause function that

$\displaystyle\mathbb{E}_{\mathbf{s}}[2^{-k}\prod_{i=1}^{n}(1-c_{m,i}s_{i})]$
$\displaystyle=0\cdot\frac{2^{n}-2^{n-k}}{2^{n}}+1\cdot\frac{2^{n-k}}{2^{n}}$
(17a) $\displaystyle=2^{-k}\,,$ (17b)

which gives the final expression

$\mathbb{E}_{\mathbf{s}}[V(\mathbf{s})]=M(1-2^{-k})\,.$ (18)

Substituting these expressions into Equation (15) gives

$\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})V(\mathbf{s})]=M(1-2^{-k})-2\frac{ML}{2^{n}}\,.$
(19)

This process needs to be iterated for the full AA effect. Fortunately, the
iterations continue to admit a closed form expression with no additional
assumptions. It can be shown that for an initial objective function
$f_{0}(\mathbf{s})$ the following relationships hold. First we define the
evolved function after each application of the pair of oracle and diffusion
operations:

$f_{\ell}(\mathbf{s})=2\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})f_{\ell-1}(\mathbf{s})]-T(\mathbf{s})f_{\ell-1}(\mathbf{s})\,.$
(20)

Letting $W_{\ell}(\mathbf{s})=T(\mathbf{s})f_{\ell}(\mathbf{s})$ we have,
equivalently,

$f_{\ell}(\mathbf{s})=2\mathbb{E}_{\mathbf{s}}[W_{\ell-1}(\mathbf{s})]-W_{\ell-1}(\mathbf{s})\,.$
(21)

The expression for $\mathbb{E}_{\mathbf{s}}[W_{0}]$ is given in Equation (15).
In general, we have

$\mathbb{E}[W_{\ell}(\mathbf{s})]=2\mathbb{E}_{\mathbf{s}}[W_{\ell-1}]\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})]-\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})W_{\ell-1}]\,,$
(22)

or, equivalently, for $\ell\geq 2$,

$\mathbb{E}_{\mathbf{s}}[W_{\ell}(\mathbf{s})]=2\mathbb{E}[W_{\ell-1}(\mathbf{s})]\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})]-\mathbb{E}_{\mathbf{s}}[W_{\ell-2}(\mathbf{s})]$
(23)

Equation (23) is obtained by expanding $\mathbb{E}[W_{\ell}(\mathbf{s})]$
using Equation (20) and applying linearity of expectation, that is,

$\begin{split}\mathbb{E}_{\mathbf{s}}\biggl{[}T(\mathbf{s})&\bigr{[}2\mathbb{E}_{\mathbf{s}}[W_{\ell-1}(\mathbf{s})]-W_{\ell-1}(\mathbf{s})\bigl{]}\biggr{]}\\\
&=2\mathbb{E}_{\mathbf{s}}[W_{\ell-1}(\mathbf{s})]\mathbb{E}[T(\mathbf{s})]-\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})W_{\ell-1}(\mathbf{s})]\,.\end{split}$
(24)

We then use the fact that

$\displaystyle\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})W_{\ell-1}(\mathbf{s})]$
$\displaystyle=\mathbb{E}_{\mathbf{s}}[T(\mathbf{s})^{2}[2\mathbb{E}_{\mathbf{s}}[W_{\ell-2}]-W_{\ell-2}]]$
(25) $\displaystyle=\mathbb{E}_{\mathbf{s}}[W_{\ell-2}(\mathbf{s})]$ (26)

since $T(\mathbf{s})^{2}=1$ at any vertex. This gives a method for computing
all of the required expectations of the form
$\mathbb{E}_{\mathbf{s}}[W_{\ell}(\mathbf{s})]$ in advance using only
knowledge of $L,n$, and $M$, and known averages. The expectations can then be
substituted as required in computations of the iterates in Equation (21).

##  IV Conclusion

(a) $f_{0}(\mathbf{s})$ at vertices  (b) Costs after 1 iteration

(c) Costs after 2 iterations  (d) Costs after 3 iterations

Figure 1:  (color online) Graphical depiction of costs for an example $k$-SAT
problem. The costs are evaluated at each of the vertices in $\mathcal{V}$,
which correspond to the costs at binary solutions, and displayed sequentially.
As such they are labeled with their binary equivalent representation. For this
example, we generated a $3$-SAT problem with $n=4$ variables, $L=1$ solution,
and $M=25$ clauses. We used the objective function formulation $V(\mathbf{s})$
to serve as $f_{0}(\mathbf{s})$, and we used the oracle function of Equation
(9) to evaluate the iterates.

Figure 1 demonstrates the transformed costs at the vertices that result from
applying conditioning to an example problem. The number of iterations that
provides the most amplification is consistent with theoretical results about
AA, i.e. $\ell^{*}=\lfloor\pi/4\sqrt{2^{n}/L}\rfloor$ iterations [13]. For
this problem, the number of iterations to apply for maximum amplification was
$3$. The most prominent effect of the conditioning is to greatly increase the
cost at $f(\mathbf{s}^{*})$ relative to the cost at all non-satisfying
vertices. The effects have discernible impact on a variety of parameters
relevant to non-convex solvers, but incur the processing of applying
AmplifySAT, for which the optimal number of iterations appears to grow as
$\mathcal{O}(\sqrt{N})$. Although, as Figure 1b and 1c highlight, there is
still a significant amplification achievable when performing fewer than the
ideal number of iterations, and the potential for a tradeoff emerges between
the cost of conditioning and the degree of amplification desired.

In simulated annealing, the maximum depth of any local minimum determines the
cooling schedule of the temperature parameter, and also impacts the
convergence rate [35, 36]. It appears that an energy landscape can be
engineered from $f_{\ell^{*}}(\mathbf{s})$ with local minima of reduced depth,
arising directly from the effects of the conditioning and potentially yielding
a faster convergence rate to the global optimum when using simulated
annealing. Similarly, solvers based on the gradient of the objective function
would use $\nabla f_{\ell^{*}}(\mathbf{s})$ after applying the algorithm, as
opposed to $\nabla f_{0}(\mathbf{s})$. A more detailed analysis of how the
gradient evolves with AmplifySAT will help illuminate whether the conditioned
gradient improves convergence properties to the global minimum. It is not
clear whether the conditioning would improve the performance of these solvers,
as the choice of oracle function $T(\mathbf{s})$ could introduce additional
local minima on the interior, or create other undesirable features across the
optimization landscape. In this sense an informed choice of oracle function
$T(\mathbf{s})$ could also greatly determine the usefulness of AmplifySAT with
gradient-based solvers.

Also of interest is whether an AmplifySAT-conditioned function can be
constructed as a continuous-time dynamical system, as was achieved with the
original UniSAT objective function [30, 22]. A successful conversion would
require another reformulation of the objective function, as the dynamical
systems are energy-minimizing (dissipative), while much of this work has
focused on amplifying, that is, increasing global maxima. However, if it is
possible to express the objective function as an energy-minimization problem,
while also improving its energy landscape through AmplifySAT, then this could
improve the performance of such solvers. The iterative nature of AmplifySAT
raises the possibility of a tradeoff between conditioning steps of the
objective function and the available energy budget.

In summary, a more comprehensive analysis will be useful to analyze the
performance of various non-convex solvers when paired with AmplifySAT. A
numerical investigation of the performance of solvers using AmplifySAT would
help delineate any performance differences and would inform theoretical
developments. While there is not currently a direct theoretical relationship
between the conditioning applied by AmplifySAT and the performance of these
solvers, it appears that the algorithm directly modifies parameters of the
objective function that affect convergence to the global optima. The
performance improvement, if it exists, will heavily depend on design decisions
regarding the original objective function, oracle function, the number of
iterations to apply, and the non-convex solution method pursued.

## Acknowledgment

The authors would like to thank The Office of Naval Research for their support
of this work under Grant No. N00014-17-1-2107, as well as Applied Research
Laboratories, The University of Texas at Austin for an internal research and
development grant.

## References

  * [1] E. Tang, “A quantum-inspired classical algorithm for recommendation systems,” in _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , 2019, pp. 217–228.
  * [2] D. Dong, C. Chen, J. Chu, and T.-J. Tarn, “Robust quantum-inspired reinforcement learning for robot navigation,” _IEEE/ASME transactions on mechatronics_ , vol. 17, no. 1, pp. 86–97, 2010.
  * [3] B. Barak, A. Moitra, R. O’Donnell, P. Raghavendra, O. Regev, D. Steurer, L. Trevisan, A. Vijayaraghavan, D. Witmer, and J. Wright, “Beating the random assignment on constraint satisfaction problems of bounded degree,” _arXiv preprint arXiv:1505.03424_ , 2015.
  * [4] N.-H. Chia, A. Gilyén, T. Li, H.-H. Lin, E. Tang, and C. Wang, “Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning,” in _Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing_ , 2020, pp. 387–400.
  * [5] A. Narayanan and M. Moore, “Quantum-inspired genetic algorithms,” in _Proceedings of IEEE international conference on evolutionary computation_. IEEE, 1996, pp. 61–66.
  * [6] N.-H. Chia, H.-H. Lin, and C. Wang, “Quantum-inspired sublinear classical algorithms for solving low-rank linear systems,” _arXiv preprint arXiv:1811.04852_ , 2018.
  * [7] Z. Chen, Y. Li, X. Sun, P. Yuan, and J. Zhang, “A quantum-inspired classical algorithm for separable non-negative matrix factorization,” _arXiv preprint arXiv:1907.05568_ , 2019.
  * [8] E. Tang, “Quantum-inspired classical algorithms for principal component analysis and supervised clustering,” _arXiv preprint arXiv:1811.00414_ , 2018\.
  * [9] L. K. Grover, “Quantum mechanics helps in searching for a needle in a haystack,” _Physical Review Letters_ , vol. 79, pp. 325–328, 1997\.
  * [10] A. Ahuja and S. Kapoor, “A quantum algorithm for finding the maximum,” _arXiv preprint quant-ph/9911082_ , 1999.
  * [11] D. S. Abrams and C. P. Williams, “Fast quantum algorithms for numerical integrals and stochastic processes,” _arXiv preprint quant-ph/9908083_ , 1999\.
  * [12] L. K. Grover, “A framework for fast quantum mechanical algorithms,” in _Proceedings of the thirtieth annual ACM symposium on Theory of computing_ , 1998, pp. 53–62.
  * [13] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation,” _Contemporary Mathematics_ , vol. 305, pp. 53–74, 2002\.
  * [14] M. W. Moskewicz, C. F. Madigan, Y. Zhao, L. Zhang, and S. Malik, “Chaff: Engineering an efficient SAT solver,” in _Proceedings of the 38th annual Design Automation Conference_ , 2001, pp. 530–535.
  * [15] N. Eén and N. Sörensson, “An extensible sat-solver,” in _International conference on theory and applications of satisfiability testing_. Springer, 2003, pp. 502–518.
  * [16] J. P. Marques-Silva and K. A. Sakallah, “Grasp: A search algorithm for propositional satisfiability,” _IEEE Transactions on Computers_ , vol. 48, no. 5, pp. 506–521, 1999.
  * [17] G. Dueck, “New optimization heuristics: The great deluge algorithm and the record-to-record travel,” _Journal of Computational physics_ , vol. 104, no. 1, pp. 86–92, 1993.
  * [18] T. Brueggemann and W. Kern, “An improved deterministic local search algorithm for 3-sat,” _Theoretical Computer Science_ , vol. 329, no. 1-3, pp. 303–313, 2004.
  * [19] S. Seitz, M. Alava, and P. Orponen, “Focused local search for random 3-satisfiability,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2005, no. 06, p. P06006, 2005.
  * [20] W. M. Spears, “Simulated annealing for hard satisfiability problems.” _Cliques, Coloring, and Satisfiability_ , vol. 26, pp. 533–558, 1993.
  * [21] U. Schöning, “A probabilistic algorithm for k-SAT and constraint satisfaction problems,” in _40th Annual Symposium on Foundations of Computer Science (Cat. No. 99CB37039)_. IEEE, 1999, pp. 410–414.
  * [22] X. Yin, B. Sedighi, M. Varga, M. Ercsey-Ravasz, Z. Toroczkai, and X. S. Hu, “Efficient analog circuits for boolean satisfiability,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 26, no. 1, pp. 155–167, 2017.
  * [23] J. D. Davis, Z. Tan, F. Yu, and L. Zhang, “Designing an efficient hardware implication accelerator for sat solving,” in _International Conference on Theory and Applications of Satisfiability Testing_. Springer, 2008, pp. 48–62.
  * [24] D. A. Basford, J. M. Smith, R. J. Connor, B. J. MacLennan, and J. Holleman, “The impact of analog computational error on an analog boolean satisfiability solver,” in _2016 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 2016, pp. 2503–2506.
  * [25] J. Gu, “Global optimization for satisfiability (SAT) problem,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 6, no. 3, pp. 361–381, 1994\.
  * [26] J. L. Johnson, “A neural network approach to the 3-satisfiability problem,” _Journal of Parallel and Distributed Computing_ , vol. 6, no. 2, pp. 435–449, 1989.
  * [27] A. Kyrillidis, A. Shrivastava, M. Vardi, and Z. Zhang, “FourierSAT: A Fourier expansion-based algebraic framework for solving hybrid boolean constraints,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 02, 2020, pp. 1552–1560.
  * [28] J. Gu, P. W. Purdom, J. Franco, and B. W. Wah, “Algorithms for the satisfiability (SAT) problem: A survey,” Cincinnati Univ oh Dept of Electrical and Computer Engineering, Tech. Rep., 1996.
  * [29] J. Gu, Q. Gu, and D. Du, “On optimizing the satisfiability (SAT) problem,” _Journal of Computer Science and Technology_ , vol. 14, no. 1, pp. 1–17, 1999\.
  * [30] M. Ercsey-Ravasz and Z. Toroczkai, “Optimization hardness as transient chaos in an analog approach to constraint satisfaction,” _Nature Physics_ , vol. 7, no. 12, pp. 966–970, 2011.
  * [31] B. R. La Cour and G. E. Ott, “Signal-based classical emulation of a universal quantum computer,” _New Journal of Physics_ , vol. 17, p. 053017, 2015.
  * [32] B. R. La Cour and C. I. Ostrove, “Subspace projection method for unstructured searches with noisy quantum oracles using a signal-based quantum emulation device,” _Quantum Information Processing_ , vol. 16, p. 7, 2017.
  * [33] Y. M. Karandashev and B. Kryzhanovsky, “Binary optimization: efficient increasing of global minimum basin of attraction,” _Optical Memory and Neural Networks_ , vol. 19, no. 2, pp. 110–125, 2010.
  * [34] I. Karandashev and B. Kryzhanovsky, “Attraction area of minima in quadratic binary optimization,” _Optical Memory and Neural Networks_ , vol. 23, no. 2, pp. 84–88, 2014.
  * [35] B. Hajek, “Cooling schedules for optimal annealing,” _Mathematics of operations research_ , vol. 13, no. 2, pp. 311–329, 1988.
  * [36] D. Mitra, F. Romeo, and A. Sangiovanni-Vincentelli, “Convergence and finite-time behavior of simulated annealing,” _Advances in applied probability_ , vol. 18, no. 3, pp. 747–771, 1986.
  * [37] R. O’Donnell, _Analysis of boolean functions_. Cambridge University Press, 2014.