revtex4-1Repair the float

††thanks: Corresponding author.

# Generating large disordered stealthy hyperuniform systems with ultra-high
accuracy to determine their physical properties

Peter K. Morse<EMAIL_ADDRESS>Department of Chemistry, Princeton
University, Princeton, NJ 08544  Department of Physics, Princeton University,
Princeton, NJ 08544  Princeton Institute of Materials, Princeton University,
Princeton, NJ 08544  Jaeuk Kim  Department of Chemistry, Princeton University,
Princeton, NJ 08544  Department of Physics, Princeton University, Princeton,
NJ 08544  Princeton Institute of Materials, Princeton University, Princeton,
NJ 08544  Paul J. Steinhardt  Department of Physics, Princeton University,
Princeton, NJ 08544  Salvatore Torquato  Department of Chemistry, Princeton
University, Princeton, NJ 08544  Department of Physics, Princeton University,
Princeton, NJ 08544  Princeton Center for Theoretical Science, Princeton
University, Princeton, NJ 08544  Princeton Institute of Materials, Princeton
University, Princeton, NJ 08544  Program in Applied and Computational
Mathematics, Princeton University, Princeton, NJ 08544

(August 28, 2024)

###### Abstract

Hyperuniform many-particle systems are characterized by a structure factor
$S({\mathbf{k}})$ that is precisely zero as
$\absolutevalue{\mathbf{k}}\rightarrow 0$; and stealthy hyperuniform systems
have $S({\mathbf{k}})=0$ for the finite range
${0<\absolutevalue{{\mathbf{k}}}\leq K}$, called the “exclusion region.”
Through a process of collective-coordinate optimization, energy-minimizing
disordered stealthy hyperuniform systems of moderate size have been made to
high accuracy, and their novel physical properties have shown great promise.
However, minimizing $S(\mathbf{k})$ in the exclusion region is computationally
intensive as the system size becomes large. In this Letter, we present an
improved methodology to generate such states using double-double precision
calculations on GPUs that reduces the deviations from zero within the
exclusion region by a factor of approximately $10^{30}$ for systems sizes more
than an order of magnitude larger. We further show that this ultra-high
accuracy is required to draw conclusions about their corresponding
characteristics, such as the nature of the associated energy landscape and the
presence or absence of Anderson localization, which might be masked, even when
deviations are relatively small.

Introduction.— Disordered hyperuniform systems are isotropic like a glass but
have long-range correlations that result in the anomalous suppression of
density fluctuations over large distances Torquato and Stillinger (2003). This
feature endows them with unique physical properties not possible in ordered
(periodic or quasiperiodic) systems Wu  _et al._ (2017); Salvalaglio  _et al._
(2020); Yu  _et al._ (2021); Yu (2023). They have been shown to be
surprisingly ubiquitous; they include perfect glasses, fermionic point
processes, disordered jammed particle packings, quantum states, certain
plasmas, galaxy distributions, eigenvalues of random matrices, and a myriad of
other examples (see Ref. Torquato (2018) and references therein).

A defining feature of a disordered hyperuniform system in $d$-dimensional
Euclidean space $\mathbb{R}^{d}$, aside from its isotropy, is that the
structure factor $S({\mathbf{k}})$ vanishes as the wavenumber $k\equiv|{{\bf
k}}|$ tends to zero Torquato and Stillinger (2003); Torquato (2018). An
important subclass is the stealthy disordered hyperuniform system in which
$S({\mathbf{k}})=0$ for a finite range of wavenumbers ${0<k\leq K}$. This
range defines the exclusion region in Fourier space where no single-scattering
events can occur Uche  _et al._ (2004, 2006); Batten  _et al._ (2008);
Torquato  _et al._ (2015), leading to them being theoretically analyzed as
$d$-dimensional “hard-sphere fluids” in Fourier space Torquato  _et al._
(2015). Disordered stealthy hyperuniform systems stand out among all non-
stealthy hyperuniform ones because they anomalously suppress density
fluctuations not only at infinite wavelengths but down to intermediate
wavelengths Torquato (2018). Moreover, “holes” in any disordered stealthy
hyperuniform point pattern are strictly bounded with a well-defined maximal
size in the thermodynamic limit Zhang  _et al._ (2017); Ghosh and Lebowitz
(2018); Torquato (2018). All of these remarkable attributes are responsible
for the novel wave, transport, and mechanical properties of disordered
stealthy hyperuniform systems Torquato (2018). For example, a scheme that
mapped stealthy point patterns onto disordered dielectric networks resulted in
the first moderately-sized amorphous photonic solid samples with complete
photonic band gaps comparable to those in periodic networks with the advantage
that the band gaps are isotropic Florescu  _et al._ (2009). Since this work,
disordered stealthy systems have been extensively fabricated Man  _et al._
(2013); Leseur  _et al._ (2016); Zhang  _et al._ (2016); Wu  _et al._ (2017);
Froufe-Pérez  _et al._ (2017); Chen and Torquato (2018); Gorsky  _et al._
(2019); Sheremet  _et al._ (2020); Kim and Torquato (2020); Tavakoli  _et al._
(2022); Gkantzounis  _et al._ (2017); Romero-García  _et al._ (2019);
Rohfritsch  _et al._ (2020); Romero-García  _et al._ (2021); Yu  _et al._
(2021); Chéron  _et al._ (2022); Klatt  _et al._ (2022) and studied
computationally especially due to their novel optical Leseur  _et al._ (2016);
Wu  _et al._ (2017); Chen and Torquato (2018); Sheremet  _et al._ (2020); Kim
and Torquato (2020); Tavakoli  _et al._ (2022), photonic Man  _et al._ (2013);
Froufe-Pérez  _et al._ (2017); Gorsky  _et al._ (2019); Yu  _et al._ (2021);
Klatt  _et al._ (2022), electronic Torquato and Chen (2018); Granchi  _et al._
(2022), phononic Gkantzounis  _et al._ (2017); Romero-García  _et al._ (2019);
Rohfritsch  _et al._ (2020); Romero-García  _et al._ (2021), and transport
Zhang  _et al._ (2016); Chéron  _et al._ (2022) properties.

Since hyperuniformity is a large-scale property of a system, it is imperative
to be able to generate sample sizes much larger than those presently possible
Zhang and Torquato (2020) in order to ascertain whether their novel properties
persist as the sample size becomes large. Indeed, this question was recently
explored in the context of photonics, where numerical evidence was presented
to show that among various disordered nonhyperuniform and hyperuniform network
dielectric solids, only certain stealthy hyperuniform ones may form bandgaps
in the thermodynamic limit Klatt  _et al._ (2022). Importantly, the bandgaps
for networks with near stealthiness (small but positive $S({\bf k})$ within
the exclusion region) eventually closed as the system size grew with a rate
inversely proportional to the distance to stealthiness $S_{\mathrm{max}}$,
which we define in the present work as

$S_{\mathrm{max}}\equiv\max_{0<k\leq K}S(\mathbf{k}).$ (1)

The work reported in Ref. Klatt  _et al._ (2022) emphasizes the importance of
computationally creating appreciably larger disordered stealthy hyperuniform
systems, but ones with the smallest value of $S_{\mathrm{max}}$ within the
exclusion region.

It is a remarkable fact that disordered stealthy hyperuniform point patterns
are highly degenerate ground states for certain non-trivial oscillatory pair
potentials (see Ref. Torquato  _et al._ (2015) and references therein). Such
ground states have been generated using a collective-coordinate energy
optimization scheme Fan  _et al._ (1991); Uche  _et al._ (2004); Batten  _et
al._ (2008); Zhang  _et al._ (2015). As such, they inevitably have a small
positive residual $S(\mathbf{k})$ in the exclusion region, i.e., small
$S_{\mathrm{max}}$, due to the level of precision with which they were
prepared. In spatial dimensions $d=1$-$3$, the largest simulated disordered
stealthy hyperuniform systems have had $N=10^{3}$, $10^{4}$, and $8\times
10^{3}$ particles, respectively Torquato  _et al._ (2021), and the highest
accuracy achieved has been $S_{\mathrm{max}}\approx 10^{-22}$ Zhang  _et al._
(2015). In order to extrapolate reliably to the thermodynamic limit,
significantly larger systems with a small $S_{\mathrm{max}}$ are required.
However, the computational cost of the standard collective-coordinate
minimization schemes scales as at least $\mathcal{O}(N^{2})$ (as derived
below), making it nearly impossible to reach the requirements without
fundamentally modifying the techniques used.

In this Letter, we show that the collective-coordinate minimization is highly
parallelizable and requires very little memory access, making it ideally
suited for a Graphical Processing Unit (GPU) based algorithm. The high degree
of parallelization afforded by GPUs allows us to create systems $20$ times
larger than previous best efforts Torquato  _et al._ (2021) in comparable time
and with $S_{\mathrm{max}}$ reduced by $30$ orders of magnitude, i.e.
$S_{\mathrm{max}}\approx 10^{-51}$ in dimensions $d=1$-$4$. This enables the
exploration of characteristics of the energy landscape associated with the
stealthy hyperuniform potential and lays the groundwork for extrapolations to
both the thermodynamic and high-dimensional ($d\to\infty$) limit. While the
techniques used are independent of dimension, in this Letter we focus on $d=2$
with system size $N=2\times 10^{6}$ as a trial case for its ease of
representation, and we use random initial conditions—corresponding to a high
temperature—to show the generic character of these results. These systems can
then be compared to the low-temperature theory of stealthy systems Torquato
_et al._ (2015) via the $\tau$ order metric, which characterizes their degree
of translational disorder.

Moreover, we demonstrate the importance of creating states with the smallest
distance to stealthiness by showing that values of $S_{\mathrm{max}}$ which
were deemed relatively small in the past can still vastly degrade certain
desired physical properties of a stealthy hyperuniform material, in addition
to those reported in Ref. Klatt  _et al._ (2022). For each application, it is
then necessary to ask what level of $S_{\mathrm{max}}$ is sufficient. In
particular, we show that for two-dimensional two-phase media derived from
stealthy hyperuniform point patterns, a full transparency interval for light
propagation in which there is no Anderson localization depends on having
$S_{\mathrm{max}}$ precisely zero over the excluded region. Thus being able to
generate stealthy hyperuniform systems over an exponentially wide range of
$S_{\mathrm{max}}$ is essential for extrapolating to the ideal condition in
the thermodynamic limit. This degree of accuracy is also necessary in
determining key properties of the ground state manifold of the stealthy
hyperuniform potential—namely its connectedness and dimension—especially near
critical points.

Definitions. The collective-coordinate method for generating disordered
stealthy hyperuniform point patterns is based on finding the ground state for
a system of point-particles whose total potential energy is the sum of pair-
wise potentials $v({\bf r})$, where $v({\bf r})$ is bounded and integrable
such that its Fourier transform $\tilde{v}({\bf k})$ exists and is a positive
function with compact support over the interval $0<|{\bf k}|\leq K$. Given a
set of $N$ points at positions $\mathbf{r}_{j}$ within a periodic box $F$ of
volume $v_{F}$ in $\mathbb{R}^{d}$, The total potential energy
$\Phi(\mathbf{r}^{N})$ has the Fourier representation Torquato  _et al._
(2015):

$\Phi(\mathbf{r}^{N})=\frac{N}{2v_{F}}\quantity[\sum_{0<k\leq
K}\tilde{v}(\mathbf{k})\mathinner{\mathcal{S}\mathopen{\left({\bf
k}\right)}}-\sum_{0<k\leq K}\tilde{v}(\mathbf{k})],$ (2)

where $\mathcal{S}({\bf k})$ is the structure factor of a _single_
configuration defined by

$\mathcal{S}(\mathbf{k})=\frac{|\tilde{n}(\mathbf{k})|^{2}}{N},$ (3)

where ${\bf k}$ is a non-zero reciprocal lattice vector of $F$, and
$\tilde{n}(\mathbf{k})$ is the complex collective density variable given by

$\tilde{n}(\mathbf{k})=\sum_{j=1}^{N}\exp(-i\mathbf{k}\cdot\mathbf{r}_{j}).$
(4)

Note that the second term in Eq. (2) is independent of structure so that one
can drop this term in practice. The potential $\Phi$ is then bounded from
below by $\Phi=0$, which defines the ground state manifold within which all
states are stealthy hyperuniform. Because $\tilde{v}(\mathbf{k})$ is strictly
positive, it does not modify the definition of the ground state manifold,
though it does provide a weight function that will funnel any dynamic or
thermal processes on the $\Phi$ landscape towards specific ground states.

Figure 1: Stealthy hyperuniform point patterns of $N=2\times 10^{6}$ particles
with $\chi=0.1$, $0.2$, $0.3$, and $0.4$ at the same number density
$\rho\equiv N/v_{F}$. Each figure is only showing $1/16^{\mathrm{th}}$ of all
the data so as to allow for better visualizations. Full configuration data for
5 iterations of each $\chi$ value are included in the Supplemental Material
sup (a).

Because the periodic hypercubic box has a finite side length $L$, all
$k$-vectors are integer linear combinations of the minimum $k$-vectors with
length $k_{\mathrm{min}}=2\pi/L$, so one can count the number for which
$0<k\leq K$. Because the structure factor has inversion symmetry, i.e.
$S(\mathbf{k})=S(-\mathbf{k})$, if there were $(2M+1)$ vectors for which
$0<k\leq K$, then only $M$ are independent. We thus state that there are $M$
independent constraints on the system. Given that there are a total of
$d(N-1)$ degrees of freedom in a $d$-dimensional point pattern, the fraction
of degrees of freedom that are constrained is

$\chi=\frac{M}{d(N-1)}.$ (5)

Because $M$ is a monotonically increasing function of $K$, $\chi$ additionally
acts as a measure of the size of the exclusion region. It has been shown that
ground states with $\chi\geq 1/2$ in dimensions for $d\leq 4$ always form
ordered hyperuniform structures, so we restrict ourselves in this study of
$d=2$ patterns to $\chi<1/2$, for which the ordered structures within the
ground state manifold are of measure zero Torquato  _et al._ (2015). Results
for one-dimensional as well as three- and higher-dimensional cases will be
reported in a future study.

The energy landscape for the stealthy hyperuniform potential in the disordered
regime ($\chi<1/2$ or $\chi<1/3$ in $d=1$) is remarkably simple. Previous
studies have shown that the inherent-structure state for any initial
condition—defined as the nearest local energy minimum accessible without a
barrier crossing Stillinger and Weber (1984)—is itself one of the many
degenerate global minima Uche  _et al._ (2004, 2006); Batten  _et al._ (2008);
Zhang  _et al._ (2015). Thus, a simple quench, achieved numerically via
standard local energy minimization techniques, will generically produce a
ground state. While previous studies found this to be true to an assigned (but
arbitrary) threshold of $\Phi<10^{-22}$ Zhang  _et al._ (2015), this
assumption will be verified to greater accuracy in this work.

Simulations. Generating disordered stealthy hyperuniform systems reduces to
finding true ground states of the total potential energy $\Phi$ in Eq. (2).
When $\chi<1/2$, the spectrum of disordered ground states is highly degenerate
with an extent that depends on the value of $\chi$, becoming infinitely
degenerate in the thermodynamic limit Torquato  _et al._ (2015). In the
collective-coordinate optimization procedure, both the initial conditions and
the choice of $\tilde{v}(\mathbf{k})$ affect the particular ground state
obtained Zhang  _et al._ (2015), though the methods used are identical. For
the purposes of illustration, we use random initial conditions defined by a
Poisson point process and the simplest form of $\tilde{v}(\mathbf{k})$. Here
$\tilde{v}(\mathbf{k})=\varepsilon_{0}L^{d}\Theta(K-k)$, where $\Theta(x)$ is
the Heaviside function and $\varepsilon_{0}$ sets the units of energy.

In order to minimize Eq. (2), we need its value and the gradient

$\mathbf{F}_{j}=-\nabla_{j}\Phi(\mathbf{r}^{N})=\frac{1}{v_{F}}\sum_{0<k\leq
K}\mathbf{k}\tilde{v}(\mathbf{k})\mathrm{Im}[\tilde{n}(\mathbf{k})\exp(i\mathbf{k}\cdot\mathbf{r}_{j})],$
(6)

which can be used as inputs into a minimizer. Here we use FIRE minimization
Bitzek  _et al._ (2006), because it requires fewer computations and converges
faster than other comparable minimizers.

By far, the most computationally expensive process of the minimization is the
calculation of the gradient, which contains two sums. First, it is necessary
to calculate all $M$ independent values of $\tilde{n}(\mathbf{k})$, each of
which contains $N$ terms. Then, one must calculate the gradient for all $N$
particles, each of which contains a sum over $M$ terms. There are thus 2
calculations—each with ${MN=\chi N(N-1)d}$ terms—which need to be recalculated
at every step of the minimization, setting the minimum time scale for
minimization as $\mathcal{O}(N^{2})$. The substantial numerical speedup comes
from the realization that both Eq. (4) and Eq. (6) are simple parallel
reduction sums, which are ideally suited for GPUs. To implement these
quenches, we use the GPU-based packing software pyCudaPack Charbonneau  _et
al._ (2012); Morse and Corwin (2014); Charbonneau  _et al._ (2015); Morse and
Corwin (2016, 2017); Arceri and Corwin (2020); Morse  _et al._ (2021), because
of its highly-optimized modular implementation of a variety of pair
interactions Charbonneau  _et al._ (2012); Morse and Corwin (2014); Arceri and
Corwin (2020); Morse  _et al._ (2021) using various minimizers Morse and
Corwin (2014) with tunable precision using double-doubles Charbonneau  _et
al._ (2015); Morse and Corwin (2016, 2017). Minimizations proceed until $\Phi$
is at the minimum value attainable within machine precision, which for
$N=2\times 10^{6}$ in unitless form is ${\Phi/\varepsilon_{0}\approx
10^{-46}}$.

It is well-established that short- and intermediate-range order increase with
the size of the exclusion region between $\chi=0$ and $1/2$ in $d=2$ Uche  _et
al._ (2004, 2006); Batten  _et al._ (2008); Zhang  _et al._ (2015); Torquato
_et al._ (2015). A positive order metric that measures the degree to which
translational order increases with $\chi$ across length scales is Torquato
_et al._ (2015)

$\tau\equiv\frac{1}{D^{d}}\int_{\mathbb{R}^{d}}\mathinner{h^{2}\mathopen{\left(r\right)}}\differential{{\bf
r}}=\frac{1}{(2\pi
D)^{d}\rho^{2}}\int_{\mathbb{R}^{d}}[S(\mathbf{k})-1]^{2}\differential{{\bf
k}},$ (7)

where $\rho$ is the number density in the thermodynamic limit,
$\mathinner{h\mathopen{\left(r\right)}}$ is the total correlation function tot
, $\mathinner{S\mathopen{\left({\bf k}\right)}}$ is the ensemble average of
Eq. (3) in the thermodynamic limit, and $D$ is a characteristic length scale,
which we take to be $D=K^{-1}$. For an ideal gas (spatially uncorrelated
Poisson point process), $\tau=0$ because
$\mathinner{h\mathopen{\left(r\right)}}=0$ for all $r$. Thus, a deviation of
$\tau$ from zero measures translational order with respect to the fully
uncorrelated case. While $\tau$ diverges for any perfect crystal and
quasicrystal in the infinite-system limit, its rate of growth as a function of
system size still provides a useful measure of the translational order in such
ordered systems Torquato  _et al._ (2018, 2019).

Figure 2: (a) The distance to stealthiness is $S_{\mathrm{max}}\approx
10^{-51}$ as stated, which is below the threshold for what can be
distinguished from zero to within machine precision. (b) The order metric
$\tau$, defined by (7), plotted as a function of $\chi$ for the simulation
data (red circles) compared to the theoretical curve from the theory of
entropically favored states Torquato  _et al._ (2015).

Results and Discussion.— We create five disordered stealthy hyperuniform
systems in $d=2$ at each value of $\chi=0.1$, $0.2$, $0.3$, and $0.4$ in order
to demonstrate the ultra-high accuracy of our GPU-based algorithm. Our results
have both high-precision—in that they are done using double-doubles— and high
accuracy, as measured by exceptionally small values of $S_{\mathrm{max}}$.
Figure 1 shows one system at each value of $\chi$, though for easier
visualization, only $1/16^{\mathrm{th}}$ of the data is shown. Full
configuration data for all 20 systems are given in the Supplemental Material
sup (a). To demonstrate the level of accuracy, we plot the structure factor on
a logarithmic scale (Fig. 2(a)). The $k$-dependence of $S(0<k\leq K)$, and
$S_{\mathrm{max}}$ in particular, depend explicitly on the form of
$\tilde{v}(\mathbf{k})$, with our particular choice yielding
$S_{\mathrm{max}}=S(k_{\mathrm{min}})$.

While $S(k>K)$ has been shown to depend strongly on initial temperature (see
Ref. Zhang  _et al._ (2015) and the Supplemental Material sup (b)), in Fig.
2(b) we find that the theoretical curve for $\tau$ derived for low-temperature
initial states departs by relatively small amounts compared to the data
obtained from high-temperature states (i.e. random initial conditions). These
results imply that while the respective pointwise behaviors of the pair
correlation function and $S(\mathbf{k})$ may sometimes differ for certain
small ranges of their arguments, the integrated measure of order across length
scales for stealthy hyperuniform states from low-$T$ and high-$T$ initial
states, as measured by $\tau$, are essentially the same. Whether this remains
true in other dimensions is a subject for future work.

Figure 3: Two-dimensional media consisting of packings of identical disks of
dielectric constant $\varepsilon_{2}=1$ in a matrix of dielectric constant
$\epsilon_{1}=11.6$. The nonoverlapping disks have radius $a$, and their
packing fraction is $\phi_{2}=0.112$. The disk centers have
$\mathinner{S\mathopen{\left(k\right)}}=S_{\mathrm{max}}$ in the exclusion
region ($k\leq K)$. (a) Semi-log plot of the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$ as a function of
a dimensionless wavenumber $k/K$. (b) Semi-log plot of the imaginary part of
the effective dynamic dielectric constant
$\imaginary[\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}]$
for transverse electric (TE) polarization as a function of the dimensionless
incident wavenumber $k_{1}/K$ in the matrix phase $1$.

We now vividly demonstrate importance of the small distance to stealthiness in
stealthy hyperuniform configurations by mapping onto two-phase dielectric
media and quantifying the degree of the attenuation of electromagnetic wave
propagation through them. Specifically, we map 2D stealthy hyperuniform point
patterns that are exactly stealthy (i.e., $S_{\mathrm{max}}=0$) in the
thermodynamic limit via theoretical methods described in Ref. Torquato  _et
al._ (2015) into a distribution of disks (phase 2) of dielectric constant
$\varepsilon_{2}=1$ in a matrix (phase 1) of dielectric constant
$\varepsilon_{1}=11.6$ by circumscribing each point with identical disks of
radius $a$ without overlap Torquato (2016a); Zhang  _et al._ (2016). By a
similar analysis, we allow $S_{\mathrm{max}}$ to be a free parameter. The
resulting area fraction covered by disks is $\phi_{2}=0.112$. The _spectral
density_ $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$ of these
disks is computed from the formula Torquato (2016a),
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}=4\pi\phi_{2}[J_{1}(ka)/k]^{2}\mathinner{S\mathopen{\left(k\right)}}$,
where $J_{1}(x)$ is the Bessel function of the first kind of order 1. The
resulting $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$ is
valid in the thermodynamic limit and inherits the distance to stealthiness
$S_{\mathrm{max}}$ of the original point patterns, in which $S_{\mathrm{max}}$
is either exactly zero or a specified positive value Sk ; see Fig. 3(a). We
obtain the effective dynamic dielectric constant
$\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}$ of
such two-phase media for incident light of transverse electric (TE)
polarization and wavenumber $k_{1}$ by using a nonlocal strong-contrast
approximation (see Eq. (73) for $d=2$ in Torquato and Kim (2021)) that depends
on $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$, which
accurately captures multiple scattering effects to all orders beyond the
quasistatic regime (i.e., $0\leq k_{1}\xi\lesssim 1$, where $\xi$ is a
characteristic inhomogeneity length scale).

The key property of interest in this system is the imaginary part of
$\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}$ which
measures the degree to which the media effectively attenuate light, as shown
in Fig. 3(b). Importantly, perfect stealthy hyperuniform media with
$S_{\mathrm{max}}=0$ exhibit full transparency (i.e.,
$\imaginary[\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}]=0$)
up to a finite wavenumber in the thermodynamic limit, implying the absence of
Anderson localization McGurn  _et al._ (1993); Aegerter and Maret (2009);
Sgrignuoli  _et al._ (2022). By contrast, a stealthy hyperuniform medium with
a moderate distance to stealthiness $(0\ll S_{\mathrm{max}}<1)$ has a positive
$\imaginary[\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}]$
proportional to $S_{\mathrm{max}}$ up to $k_{1}=\mathcal{O}(K)$, implying that
transmittance through this medium is increasingly suppressed with a larger
sample size or a larger $S_{\mathrm{max}}$, and that Anderson localization
likely emerges. Therefore, to better understand the insulator-conductor
transition as the system size increases, it is essential to generate much
larger sample sizes with a small $S_{\mathrm{max}}$ via the techniques
described in this paper.

_Conclusions_ — Through the use of GPU minimizations with double-double
precision, we have been able to dramatically increase the size and reduce the
distance to stealthiness of stealthy hyperuniform point patterns. The ability
to create disordered stealthy hyperuniform systems of both large sizes and
ultra-small $S_{\mathrm{max}}$ is imperative to study their novel properties.
While the requirement for large sample sizes has been recognized, the need for
a small $S_{\mathrm{max}}$ has gone largely unnoticed due to the inability to
create large samples. Indeed, most previous studies were based on relatively
small sample sizes ($100\lesssim N\lesssim 1000$) with $10^{-10}\lesssim
S_{\mathrm{max}}\lesssim 10^{-5}$. Using a capacity to create much larger
systems with very small values of $S_{\mathrm{max}}$, we can shed light on the
emerging novel properties of stealthy hyperuniform systems that depend on
having a small distance to stealthiness. One such open question is whether
isotropic photonic (or phononic) band gaps exist in two-phase systems derived
from stealthy hyperuniform point patterns in the thermodynamic limit Klatt
_et al._ (2022). In addition, while disordered stealthy hyperuniform systems
of moderate sizes are known to be transparent up to a finite $k$ for
electromagnetic Leseur  _et al._ (2016); Froufe-Pérez  _et al._ (2017);
Torquato and Kim (2021) and elastic waves Rohfritsch  _et al._ (2020); Romero-
García  _et al._ (2021); Kim and Torquato (2020), it is not known whether this
transparency interval persists in the thermodynamic limit without Anderson
localization. Our theoretical analysis [summarized in Fig. 3(b)] has shown
that this transparency property can increasingly degrade as $S_{\mathrm{max}}$
and the system size are made larger, stressing the importance of making
$S_{\mathrm{max}}$ for a given system as small as possible, which is now
achieved with ultra-high accuracy via our improved numerical methodology.
Samples of much larger size with vastly smaller $S_{\mathrm{max}}$, heretofore
not possible, can now be potentially fabricated by combining our designs with
photolithographic and 3D printing techniques Tumbleston  _et al._ (2015);
Shirazi  _et al._ (2015); Zhao and Mason (2018) in order to explore the
presence or absence of Anderson localization McGurn  _et al._ (1993); Aegerter
and Maret (2009); Wiersma (2013) in disordered stealthy hyperuniform systems
as function of sample size.

These new high precision results also lay the groundwork to further
investigate the characteristics of the energy landscape of the stealthy
hyperuniform potential. Indeed, we now have abundant evidence that local
minima do not exist (or are extremely rare) for $\chi<1/2$, implying that
states with high values of $S_{\mathrm{max}}$ found in previous studies are
not locally stable; rather, they are the result of incomplete minimization,
making them excited states with an effective temperature and a modified
isothermal compressibility Torquato  _et al._ (2015). Future work will aim to
understand the nature of the true ground state manifold, particularly its
connectivity, the topology near crystalline states, and the topological change
that occurs when $\chi=1/2$ (especially as a function of dimension $d$) and
how these landscape properties modify the novel physical properties of
disordered stealthy hyperuniform systems. Moreover, with significantly larger
systems, one will be able to probe higher-order correlation functions beyond
$\tau$ in determining the structural order in stealthy hyperuniform systems
Torquato  _et al._ (2021). Additionally, while we have focused here on
isotropic exclusion regions, it is straightforward to apply our ultra-high
accuracy method to create larger stealthy hyperuniform systems with
anisotropic exclusion regions, and thus statistically anisotropic ground
states Martis  _et al._ (2013); Torquato (2016b).

###### Acknowledgements.

The Research was sponsored by the Army Research Office and was accomplished
under Cooperative Agreement Number W911NF-22-2-0103. Simulations were
performed on computational resources managed and supported by the Princeton
Institute for Computational Science and Engineering (PICSciE). P.K.M. would
like to thank Eric Corwin for discussions related to pyCudaPack.

## References

  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
  * Wu  _et al._ (2017) B.-Y. Wu, X.-Q. Sheng, and Y. Hao, PLoS One 12, e0185921 (2017).
  * Salvalaglio  _et al._ (2020) M. Salvalaglio, M. Bouabdellaoui, M. Bollani, A. Benali, L. Favre, J.-B. Claude, J. Wenger, P. de Anna, F. Intonti, A. Voigt, and M. Abbarchi, Phys. Rev. Lett. 125, 126101 (2020).
  * Yu  _et al._ (2021) S. Yu, C.-W. Qiu, Y. Chong, S. Torquato, and N. Park, Nat. Rev. Mater. 6, 226 (2021).
  * Yu (2023) S. Yu, Nat. Comput. Sci. 3, 128 (2023).
  * Torquato (2018) S. Torquato, Phys. Rep. 745, 1 (2018).
  * Uche  _et al._ (2004) O. U. Uche, F. H. Stillinger, and S. Torquato, Phys. Rev. E 70, 046122 (2004).
  * Uche  _et al._ (2006) O. U. Uche, S. Torquato, and F. H. Stillinger, Phys. Rev. E 74, 031104 (2006).
  * Batten  _et al._ (2008) R. D. Batten, F. H. Stillinger, and S. Torquato, J. Appl. Phys. 104, 033504 (2008).
  * Torquato  _et al._ (2015) S. Torquato, G. Zhang, and F. H. Stillinger, Phys. Rev. X 5, 021020 (2015).
  * Zhang  _et al._ (2017) G. Zhang, F. H. Stillinger, and S. Torquato, Soft Matter 13, 6197 (2017).
  * Ghosh and Lebowitz (2018) S. Ghosh and J. L. Lebowitz, Commun. Math. Phys. 363, 97 (2018).
  * Florescu  _et al._ (2009) M. Florescu, S. Torquato, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 106, 20658 (2009).
  * Man  _et al._ (2013) W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato, P. M. Chaikin, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 110, 15886 (2013).
  * Leseur  _et al._ (2016) O. Leseur, R. Pierrat, and R. Carminati, Optica 3, 763 (2016).
  * Zhang  _et al._ (2016) G. Zhang, F. H. Stillinger, and S. Torquato, J. Chem. Phys. 145, 244109 (2016).
  * Froufe-Pérez  _et al._ (2017) L. S. Froufe-Pérez, M. Engel, J. J. Sáenz, and F. Scheffold, Proc. Natl. Acad. Sci. U.S.A. 114, 9570 (2017).
  * Chen and Torquato (2018) D. Chen and S. Torquato, Acta Mater. 142, 152 (2018).
  * Gorsky  _et al._ (2019) S. Gorsky, W. A. Britton, Y. Chen, J. Montaner, A. Lenef, M. Raukas, and L. Dal Negro, APL Photonics 4, 110801 (2019).
  * Sheremet  _et al._ (2020) A. Sheremet, R. Pierrat, and R. Carminati, Phys. Rev. A 101, 053829 (2020).
  * Kim and Torquato (2020) J. Kim and S. Torquato, Proc. Natl. Acad. Sci. U.S.A. 117, 8764 (2020).
  * Tavakoli  _et al._ (2022) N. Tavakoli, R. Spalding, A. Lambertz, P. Koppejan, G. Gkantzounis, C. Wan, R. Röhrich, E. Kontoleta, A. F. Koenderink, R. Sapienza, M. Florescu, and E. Alarcon-Llado, ACS Photonics 9, 1206 (2022).
  * Gkantzounis  _et al._ (2017) G. Gkantzounis, T. Amoah, and M. Florescu, Phys. Rev. B 95, 094120 (2017).
  * Romero-García  _et al._ (2019) V. Romero-García, N. Lamothe, G. Theocharis, O. Richoux, and L. García-Raffi, Phys. Rev. Appl. 11, 054076 (2019).
  * Rohfritsch  _et al._ (2020) A. Rohfritsch, J.-M. Conoir, T. Valier-Brasier, and R. Marchiano, Phys. Rev. E 102, 053001 (2020).
  * Romero-García  _et al._ (2021) V. Romero-García, É. Chéron, S. Kuznetsova, J.-P. Groby, S. Félix, V. Pagneux, and L. García-Raffi, APL Mater. 9, 101101 (2021).
  * Chéron  _et al._ (2022) É. Chéron, S. Félix, J.-P. Groby, V. Pagneux, and V. Romero-García, Appl. Phys. Lett. 121, 061702 (2022).
  * Klatt  _et al._ (2022) M. A. Klatt, P. J. Steinhardt, and S. Torquato, Proc. Natl. Acad. Sci. U.S.A. 119, e2213633119 (2022).
  * Torquato and Chen (2018) S. Torquato and D. Chen, Multifunct. Mater. 1, 015001 (2018).
  * Granchi  _et al._ (2022) N. Granchi, R. Spalding, M. Lodde, M. Petruzzella, F. W. Otten, A. Fiore, F. Intonti, R. Sapienza, M. Florescu, and M. Gurioli, Adv. Opt. Mater. 10, 2102565 (2022).
  * Zhang and Torquato (2020) G. Zhang and S. Torquato, Phys. Rev. E 101, 032124 (2020).
  * Fan  _et al._ (1991) Y. Fan, J. K. Percus, D. K. Stillinger, and F. H. Stillinger, Phys. Rev. A 44, 2394 (1991).
  * Zhang  _et al._ (2015) G. Zhang, F. H. Stillinger, and S. Torquato, Phys. Rev. E 92, 022119 (2015).
  * Torquato  _et al._ (2021) S. Torquato, J. Kim, and M. A. Klatt, Phys. Rev. X 11, 021028 (2021).
  * sup (a) Supplementary data will be uploaded at the time of final publication.
  * Stillinger and Weber (1984) F. H. Stillinger and T. A. Weber, Science 225, 983 (1984).
  * Bitzek  _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Phys. Rev. Lett. 97, 170201 (2006).
  * Charbonneau  _et al._ (2012) P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 109, 205501 (2012).
  * Morse and Corwin (2014) P. K. Morse and E. I. Corwin, Phys. Rev. Lett. 112, 115701 (2014).
  * Charbonneau  _et al._ (2015) P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 114, 125504 (2015).
  * Morse and Corwin (2016) P. K. Morse and E. I. Corwin, Soft Matter 12, 1248 (2016).
  * Morse and Corwin (2017) P. K. Morse and E. I. Corwin, Phys. Rev. Lett. 119, 118003 (2017).
  * Arceri and Corwin (2020) F. Arceri and E. I. Corwin, Phys. Rev. Lett. 124, 238002 (2020).
  * Morse  _et al._ (2021) P. K. Morse, S. Roy, E. Agoritsas, E. Stanifer, E. I. Corwin, and M. L. Manning, Proc. Natl. Acad. Sci. U.S.A. 118, e2019909118 (2021).
  * (45) The total correlation function $h(r)$ is related to the pair correlation function $g_{2}(r)$ as $h(r)\equiv g_{2}(r)-1$.
  * Torquato  _et al._ (2018) S. Torquato, G. Zhang, and M. de Courcy-Ireland, J. Stat. Mech. 2018, 093401 (2018).
  * Torquato  _et al._ (2019) S. Torquato, G. Zhang, and M. D. Courcy-Ireland, J. Phys. A: Math. Theor. 52, 135002 (2019).
  * sup (b) See Supplementary Materials for additional plots of the pair correlation function and structure factor shown in Fig. 2.
  * Torquato (2016a) S. Torquato, J. Phys.: Condens. Matter 28, 414012 (2016a).
  * (50) Since disordered stealthy hyperuniform point patterns can be treated as hard-sphere fluids in Fourier space Torquato  _et al._ (2015), $S(k)$ decays to unity exponentially fast. Thus, the precise form of $S(k)$ outside the exclusion region is irrelevant for our purposes because the value of $S(k)$ within the exclusion region determines the value of $\mathrm{Im}[\mathinner{\varepsilon_{e}^{\mathrm{TE}}\mathopen{\left(k_{1}\right)}}]$ within the transparency interval. For this reason, we do not present the precise form here.
  * Torquato and Kim (2021) S. Torquato and J. Kim, Phys. Rev. X 11, 021002 (2021).
  * McGurn  _et al._ (1993) A. R. McGurn, K. T. Christensen, F. M. Mueller, and A. A. Maradudin, Phys. Rev. B 47, 13120 (1993).
  * Aegerter and Maret (2009) C. M. Aegerter and G. Maret, in  _Progress in Optics_ , Vol. 52 (Elsevier, 2009) pp. 1–62.
  * Sgrignuoli  _et al._ (2022) F. Sgrignuoli, S. Torquato, and L. Dal Negro, Phys. Rev. B 105, 064204 (2022).
  * Tumbleston  _et al._ (2015) J. R. Tumbleston, D. Shirvanyants, N. Ermoshkin, R. Janusziewicz, A. R. Johnson, D. Kelly, K. Chen, R. Pinschmidt, J. P. Rolland, A. Ermoshkin, E. T. Samulski, and J. M. DeSimone, Science 347, 1349 (2015).
  * Shirazi  _et al._ (2015) S. F. S. Shirazi, S. Gharehkhani, M. Mehrali, H. Yarmand, H. S. C. Metselaar, N. Adib Kadri, and N. A. A. Osman, Sci. Technol. Adv. Mater. 16, 033502 (2015).
  * Zhao and Mason (2018) K. Zhao and T. G. Mason, Rep. Prog. Phys. 81, 126601 (2018).
  * Wiersma (2013) D. S. Wiersma, Nat. Photon. 7, 188 (2013).
  * Martis  _et al._ (2013) S. Martis, É. Marcotte, F. H. Stillinger, and S. Torquato, J Stat Phys 150, 414 (2013).
  * Torquato (2016b) S. Torquato, Phys. Rev. E 94, 022122 (2016b).


# Machine Learning for SAT:
Restricted Heuristics and New Graph Representations

Mikhail Shirokikh  St. Petersburg State University, St. Petersburg, Russia
Ilya Shenbin  Steklov Institute of Mathematics at St. Petersburg,
St. Petersburg, Russia
Anton Alekseev  Steklov Institute of Mathematics at St. Petersburg,
St. Petersburg, Russia  St. Petersburg State University, St. Petersburg,
Russia  Sergey Nikolenko  Steklov Institute of Mathematics at St. Petersburg,
St. Petersburg, Russia  St. Petersburg State University, St. Petersburg,
Russia

###### Abstract

Boolean satisfiability (SAT) is a fundamental NP-complete problem with many
applications, including automated planning and scheduling. To solve large
instances, SAT solvers have to rely on heuristics, e.g., choosing a branching
variable in DPLL and CDCL solvers. Such heuristics can be improved with
machine learning (ML) models; they can reduce the number of steps but usually
hinder the running time because useful models are relatively large and slow.
We suggest the strategy of making a few initial steps with a trained ML model
and then releasing control to classical heuristics; this simplifies cold start
for SAT solving and can decrease both the number of steps and overall runtime,
but requires a separate decision of when to release control to the solver.
Moreover, we introduce a modification of _Graph-Q-SAT_ tailored to SAT
problems converted from other domains, e.g., open shop scheduling problems. We
validate the feasibility of our approach with random and industrial SAT
problems.

##  1 Introduction

Boolean satisfiability (SAT), i.e., deciding whether a Boolean formula in
conjunctive normal form (CNF) is satisfiable, is the archetypal NP-complete
problem, with numerous applications in computer science. There exist several
different classes of approaches to solving SAT: stochastic local search (SLS),
Davis–Putman–Logemann–Loveland solvers (DPLL), conflict-driven clause learning
(CDCL), ordered binary decision diagrams (OBDD), and others [10]. All known
methods are, naturally, exponential in the worst case, and to solve large
problems in practice they need to use various heuristics that have a crucial
impact on performance: SLS solvers choose a variable to flip, DPLL and CDCL
solvers choose branching variables and their assignments, OBDD-based solvers
choose clauses for conjunction and variables for projections. In all cases,
key decisions have to be made under uncertainty, which opens up possibilities
for using machine learning (ML) techniques to improve SAT and SMT solvers.

Various ML-based techniques have arisen to make better heuristic decisions
(see “Related work”). However, they introduce another tradeoff: ML models are
usually computationally heavy, and benefits in terms of the number of solver
iterations have to be weighted against extra costs the model incurs. Here, we
consider _Graph-Q-SAT_ [13], a model that can reduce the number of iterations
for a CDCL solver but that is based on a reinforcement learning (RL) agent
with a complex graph neural network (GNN) inside.

In this work, we propose a technique for finding this tradeoff: we suggest to
make only a few initial steps with a trained RL agent and then release control
to classical heuristics. We consider several modifications of this idea:
constant number of “heavy” steps, a separate action in the RL agent trained to
release control, and a separate head in the model to make this decision. We
introduce a novel modification with an action pool that uses the (heavy) RL
agent once to predict several actions that can be used for several steps.
Moreover, we propose a new approach tailored specifically for SAT instances
originating from other optimization problems such as open shop scheduling
(OSSP). Namely, we construct a GNN whose graph corresponds not to the SAT
formula but to the original OSSP instance; this greatly reduces the size of
the graph and improves the results in SAT solving. We validate our approaches
with an evaluation study on both random and industrial SAT instances.

The paper is organized as follows. Section 2 surveys related work on machine
learning approaches to SAT. Section 3 introduces _Graph-Q-SAT_ (Section 3.1)
and our novel modifications for it (Section 3.2). Section 4 reports the
experimental setup and datasets and specifies several novel implementation
modifications that further improve _Graph-Q-SAT_. Section 5 presents and
discusses the results of our evaluation study, and Section 6 concludes the
paper.

##  2 Related work

SAT is a classification problem; end-to-end machine learning (ML) approaches
have been proposed by Devlin and O’Sullivan [6] and with graph neural networks
(GNN) by Selsam et al. [24], a method later extended to Circuit-SAT by
Amizadeh et al. [1] and MAX-SAT by Liu et al. [17] and improved with recurrent
neural networks by Ozolins et al. [21]. However, while a satisfying assignment
can be verified, end-to-end solutions for an exact combinatorial problem via
statistical methods cannot replace a complete solver, i.e., cannot make sure
that a formula is unsatisfiable.

Instead of an end-to-end replacement, ML can be used to augment complete SAT
solvers by enhancing their _heuristics_. SAT solvers–Davis-Putnam-Logeman-
Loveland (DPLL), conflict-driven clause learning (CDCL), and stochastic local
search (SLS)—have decisions that have to be made during operation with
incomplete information, usually recursively. Such decisions can be improved
with ML. Boolean variables initialization has been improved by Wu et al. [33];
in SLS solvers, variable re-initialization is called upon every restart, so
Zhang et al. [35] train a neural network to improve SLS initial assignments.
Nejati et al. [20] learned restart strategies for CDCL solvers while Han et
al. [8] predicted glue variables for the _CaDiCaL_ solver previously proposed
by Biere et al. [3].

The most important direction here is _variable selection_ for branching in
DPLL and CDCL solvers and flipping in SLS solvers. For DPLL and CDCL, the goal
is to improve over the _Variable State Independent Decaying Sum_ (VSIDS)
heuristic originally introduced in the _Chaff_ solver by Moskewicz et al. [19]
that favors recently learned clauses based on specially introduced activity
scores; Biere and Frölich [4] survey VSIDS variations. Selsam et al. [23]
learn to predict the unsatisfiable core of a formula, using predictions to re-
rank variables for branching; Jaszczur et al. [12] use a GNN in a similar way,
while Han [9] prioritizes variables that often occur in DRAT proofs, learning
to predict this fact.

However, the main downside of ML heuristics for solving SAT is the trade-off
between reducing the number of steps and the time required to make predictions
with the help of the heavyweight modern machine learning models. Therefore,
while most works report significant reductions in the number of iterations
from the use of ML-improved heuristics, one should also track the actual
running time of the solver with and without these heuristics on equivalent
hardware, and most ML methods have limited usability in practice due to
computational inefficiency: a heavy neural network has to be used on every
branching step. _NeuroComb_ proposed by Wang et al. [31] uses GNNs only for
preprocessing, running inference once per input formula (_static_
predictions). Many works use reinforcement learning (RL) [27] to learn
heuristics based on the running time or number of iterations as a reward. In
particular, Yolcu and Póczos [34] use RL for SLS solvers to choose flip
variables, Han [8] extends _NeuroCore_ with RL, Liang et al. [15, 16] apply
multi-armed bandits as branching heuristics, Lagoudakis and Littman [14] use
TD-learning to select a branching heuristic given the current state, while
Nejati et al. [20] apply multi-armed bandits to choose a restart strategy.

_Graph-Q-SAT_ by Kurin et al. [13] utilizes Q-learning for branching
heuristics in CDCL solvers via a Markov decision process (MDP) whose states
are defined by unassigned variables and unsatisfied clauses with these
variables; we will describe this approach in detail below. We also note here
that Wang et al. [30] already introduced a similar MDP formulation, which did
not use GNNs and instead was based on _AlphaGo Zero_. Interestingly, games
such as Go and chess indeed give rise to computationally intractable problems,
e.g., the game of Go on an $n\times n$ board has been shown to be EXPTIME-
complete [22] and chess to be PSPACE-complete [26], so these two directions of
research are not as different as it might seem.

##  3 Proposed method

###  3.1 Graph-Q-SAT

Initialize the formula Choose a literal, assign to True Propagate Boolean
constraints Is there a conflict? Are all variables assigned? Find conflicting
condition, add its negation to the formula Output SAT Output UNSAT NoYesNoYes

Figure 1: CDCL solver operation: choosing an assignment is a heuristic step
that can be improved by machine learning.

In this work, we build upon the _Graph-Q-SAT_ approach by Kurin et al. [13]
that uses Q-learning for branching heuristics in CDCL solvers [25]. A CDCL
solver operates as shown in Fig. 1: on every step, it

  * •

chooses a literal (variable and its value) to assign,

  * •

then propagates Boolean constraints such as, e.g., unit clauses (clauses with
one literal) that lead to more necessary assignments.

If the solver arrives at a conflict (empty clause), it finds the conflicting
condition (a cut in the implication graph), backtracks to the assignment where
it occurred, and adds the clause representing its negation to the formula. If
the algorithm has backtracked from both possible assignments of a variable,
the formula is unsatisfiable. If all variables have been assigned with no
conflict, the solver has found a satisfying assignment.

On every iteration, a CDCL solver has to choose a literal. This part is a
heuristic, and choosing correct variables can lead to exponential speedups in
the running time and number of iterations needed for both SAT and UNSAT
instances, so this is the place where ML-based heuristics may help.

_Graph-Q-SAT_ is a reinforcement learning method that learns to choose
literals by treating the problem as a Markov decision process (MDP):

  * •

the task $\tau\sim
D(\phi,(\mathrm{un})\mathrm{SAT},n_{\mathrm{vars}},n_{\mathrm{clauses}})$ is a
SAT problem, sampled from a task family $\phi$ with known satisfiability
status and given number of variables and clauses;

  * •

each state of the MDP is defined by unassigned variables and unsatisfied
clauses with these variables;

  * •

actions available to the agent are the choice of a variable and its polarity;

  * •

a terminal state is reached when either a satisfying assignment or an unsat
core is found, and the final reward is $1$ with a discount factor $\gamma<1$
to reward solving SAT faster.

$x_{1}$$x_{2}$$x_{3}$$x_{4}$$c_{1}$$c_{2}$$[1,0]$$[1,0]$$[0,1]$$[0,1]$$[1,0]$$[1,0]$ | $x_{1}$$x_{2}$$x_{3}$$x_{4}$$c_{1}$$\rho_{e\to v}$$\phi_{v}$$c_{2}$$\phi_{e}$$\mathbf{u}$$\mathbf{e}_{11}$$\mathbf{e}_{21}$$\mathbf{e}_{31}$$\mathbf{v}_{c2}$$\mathbf{v}_{x4}$$\mathbf{v}_{c1}$
---|---
(a) | (b)
Figure 2: _Graph-Q-SAT_ : (a) the variable-clause graph for the formula
$(x_{1}\lor x_{2}\lor\lnot x_{3})\land(\lnot x_{2}\lor x_{3}\lor x_{4})$; (b)
computation of sample GNN updates for $\mathbf{v}_{c1}$ (top) and
$\mathbf{e}_{x4,c2}$ (bottom).

This is a sufficiently general formulation that has been only slightly
modified in other works that apply RL to SAT. _Graph-Q-SAT_ uses _Q-learning_
that learns to predict expected utilities for each action; for SAT, it means
that for a given formula, _Graph-Q-SAT_ learns to predict expected utilities
$Q(l)$ (running time rewards) based on the formula structure for each literal
$l$; we can choose $\arg\max_{l}Q(l)$ as the next literal to assign during
inference. See, e.g., Mnih et al. [18] for a detailed discussion of (deep)
Q-learning.

To estimate $Q(l)$, _Graph-Q-SAT_ uses a _graph neural network_ (GNN) [2],
which is a common theme in ML for SAT [24, 34]. It uses the variable-clause
graph to represent formulas, i.e., a formula with $n$ variables and $m$
clauses corresponds to a bipartite graph $G=(V,E)$ with $n$ and $m$ variables
in each part (see Fig. 2a). An edge $(x_{i},c_{j})\in E$ means that $x_{i}$
occurs in the clause $c_{j}$; an edge is labeled $[0,1]$ if $c_{j}$ contains
the positive literal $x_{i}$ and $[1,0]$ if it contains $\lnot x_{i}$.

On each step, GNN propagates information along the edges of $G$, updating the
embeddings of edges $\mathbf{e}_{ij}$, embeddings of vertices
$\mathbf{v}_{i}$, and a global context vector $\mathbf{u}$. The network is
defined by update functions $\phi_{e}$, $\phi_{v}$, $\phi_{u}$ that change the
embeddings and aggregation functions $\rho_{e\to v}$, $\rho_{e\to u}$,
$\rho_{v\to u}$ that combine the embeddings of adjacent entities.

Formally, on every iteration the GNN updates (see Fig. 2b)

$\displaystyle\mathbf{e}_{ij}:=$
$\displaystyle\phi_{e}\left(\mathbf{u},\mathbf{e}_{ij},\mathbf{v}_{i},\mathbf{v}_{j}\right)\quad\forall\
e_{ij}\in E,$ $\displaystyle\mathbf{v}_{i}:=$
$\displaystyle\phi_{v}\left(\mathbf{u},\mathbf{v}_{i},\rho_{e\to
v}\left(\left\\{\mathbf{e}_{ki}\mid e_{ki}\in
E\right\\}\right)\right)\quad\forall\ v_{i}\in V,$ $\displaystyle\mathbf{u}:=$
$\displaystyle\phi_{u}\left(\mathbf{u},\rho_{e\to
u}\left(\left\\{\mathbf{e}\mid e\in E\right\\}\right),\rho_{v\to
u}\left(\left\\{\mathbf{v}\mid v\in V\right\\}\right)\right),$

repeating until convergence or until it is stopped. Sample computations are
illustrated in the graph shown in Fig. 2b.

The functions are parameterized with neural networks, trained with
backpropagation on the learning phase. The original _Graph-Q-SAT_ uses the
_Encoder_ – _Core_ – _Decoder_ architecture that consists of three
significantly different parts. The encoder and decoder do not aggregate along
the edges; they apply linear layers with _ReLU_ activation functions to the
feature vectors $\mathbf{v}$, $\mathbf{e}$, and $\mathbf{u}$. The core is
located between them and consists of 4 GNN layers as described above.
Interestingly, while the original _Graph-Q-SAT_ code does have generic encoder
and decoder parts implemented, the version introduced by Kurin et al. [13] has
no hidden layers in the encoder and decoder.

As the underlying solver _Graph-Q-SAT_ uses _MiniSAT_ [7], an award-winning
CDCL solver that chooses the next branching literal based on their _activity
scores_ , with a variation of the _Variable State Independent Decaying Sum_
(VSIDS) heuristic; it favors recently learned clauses, increasing activity
scores for variables in recent clauses and decaying them with time. _Graph-Q-
SAT_ replaces VSIDS scores with predicted values of $Q(l)$.

###  3.2 Our modifications

We have made several modifications and improvements for the _Graph-Q-SAT_
framework. First and most importantly, we note that running a graph neural
network on every step of SAT solving is a very resource-intensive task; even
if it helps further reduce the number of steps a little, it certainly hurts
the running time, usually drastically. The most influential decisions come at
the beginning, and this is exactly when classical SAT solving heuristics such
as MiniSAT activity scores are not yet meaningful (they don’t yet have enough
statistics). Therefore, we consider two ways of limiting the time:

  1. (1)

running _Graph-Q-SAT_ for a fixed number of steps and then releasing control
to MiniSAT;

  2. (2)

learning when to release control to MiniSAT, i.e., training a separate head
that predicts this.

We also experimented with running the model once per several steps, but these
experiments did not bring improvements and are not reported below.

Solver state Clause-variable graph Q-values Gym env.GNN Action Variable
assignment AgentGym env.Minisat
---
(a)
Solver state Action pool empty? Clause-variable graph Q-values YesGNN Action
pool: Action 1 $\ldots$ Action $k$ Variable assignment AgentGym env.MinisatNo
(b)
Figure 3: Using the action pool: (a) regular _Graph-Q-SAT_ operation; (b) our
modification with the action pool. Green arrows show “fast” operations; red
arrows, “slow” operations that incur significant computational costs.

Second, we propose one more way to reduce the number of calls to the agent,
namely to utilize every run of the agent network in a more efficient way. The
graph that represents the current formula changes only slightly after each
step of the SAT solving agent. Therefore, we can assume that the list of
actions with the highest Q-values in most cases will also change only
slightly, which is confirmed by experiments. Following this motivation, we
propose to keep top $k$ actions that correspond to the highest Q-values and
perform these actions sequentially (ignoring actions that have become invalid
by the time of their execution). Such a policy is used at test time only, not
at train time. The number $k$, i.e., the number of actions per agent run, is a
constant hyperparameter, although we note as a possible extension for future
work that the number $k$ could also be adaptive. The value of $k$ should be
large enough to get sufficient runtime acceleration in SAT solving while
remaining small enough so that the environment does not change too much and
the values of the $Q$ function remain relevant. We call this set of
precomputed actions the _action pool_ ; the resulting workflow is shown in
Fig. 3b in comparison with the regular _Graph-Q-SAT_ workflow (Fig. 3a).

Third, we have experimented with _Graph-Q-SAT_ architectures. In the final
version, we:

  * •

turned on the encoder and decoder parts;

  * •

increased the number of core layers from 4 to 13;

  * •

increased the depth of core layers from 1 to 2;

  * •

added a separate “release control” action predicted with two linear layers
over the penultimate layer of the network.

Third, Kurin et al. [13] use a MiniSAT implementation that serves as the RL
environment for _Graph-Q-SAT_ , and this significantly increases the runtime.
Since we release control to MiniSAT and do not return to RL, we are able to
use the original MiniSAT after release, which has sped up SAT solving by a
factor of about 10x; in all experiments, we report the runtimes of our
version.

##  4 Experimental setup

###  4.1 Model setup and training

Similar to Kurin et al. [13], we used Random 3-SAT instances during training,
with $50$ variables and $218$ clauses each (close to the phase transition
threshold), and coloring problems for random graphs with $75$ vertices and
$180$ edges; see also the SATLIB benchmark [11].

Since we need to learn to switch between RL and VSIDS heuristics, we have
changed the rewards. With constant $-p$ reward for every step, as in the
original Graph-Q-SAT [13], the agent would have no motivation to flip the
switch as long as his variable assignments are better than random. Therefore,
instead of a constant penalty $-p$ for every iteration we used $-t$, where $t$
is the CPU time in milliseconds for this step (higher for RL and lower for
VSIDS).

We also turned off discounting for $Q$ values (set $\gamma=1$) to motivate the
agent to look ahead to long-term costs and switched to Double DQN instead of
regular DQN; both of these changes significantly improved the agent’s
behaviour in our experiments.

We have also expanded the action space with a new one that releases control to
VSIDS. During training, we force the agent to make a small constant number of
steps $n$ since otherwise for an untrained agent this action would always be
preferable; the value of $n$ did not have a significant effect on the results.
We also replaced DQN learning with Double DQN [29], which stabilized the
training process.

Note that the reward received after the switch varies over time, which does
not match the conventional definition of MDPs. However, we can treat this
model as an MDP as follows: suppose that at time $t$ we have a valid action
with reward $r_{t}$ that switches control to VSIDS, and at time $t+1$ we have
a different action that does the same but obtains reward $r_{t+1}$. In the
implementation, it turns out that parameters of the action with reward
$r_{t+1}$ are initialized by parameters of the action with reward $r_{t}$.

Each model is run three times on every task, and the tables show averaged
results.

###  4.2 Evaluation datasets

For evaluation, we have chosen several traditional SAT datasets and generated
open shop scheduling problems (OSSP).

####  $\mathbf{SR}(n)$.

For random generation, we use the $\mathbf{SR}(n)$ distribution introduced by
Selsam et al. [24]. To generate a random clause on $n$ variables,
$\mathbf{SR}(n)$ samples a small integer $k$ with mean a little over $4$,
namely as a sum of two random variables

$1+\mathbf{Bernoulli}(0.7)+\mathbf{Geo}(0.4)$

(since there should not be many $2$-literal clauses), then samples $k$
variables uniformly at random without replacement, and negates each one with
probability $0.5$.

The clauses $c_{i}$ are consecutively generated and added to the SAT instance
until the problem becomes unsatisfiable at step $m$. Since we know that it is
only $c_{m}$ that “breaks” satisfiability, negating any literal in $c_{m}$
leads to a satisfiable formula with clauses
$\\{c_{1},...,c_{m-1},c^{\prime}_{m}\\}$; both $c_{1}\land...\land
c_{m-1}\land c^{\prime}_{m}$ and $c_{1}\land...\land c_{m}$ are added to the
set. This is a standard distribution for SAT evaluation [24, 34].

#### Industrial.

This is the _maris05_ subset of industrial problems from the SAT 2005
Competition benchmark111Available at https://zenodo.org/record/6528885.

#### OSSP datasets.

Open shop scheduling is a standard scheduling problem; the task is to
distribute the workload across a finite number of resources (machines). The
workload is represented as a list of jobs; a job is essentially a set of
order-independent operations that each must be processed on a certain machine,
i.e., we need to schedule $n$ jobs on $m$ machines, where each of the $n$ jobs
has its own sequence of operations that have to be done within $m$ machines,
with the goal to optimize an objective function under constraints that can
define termination times, delay times, or total flow times, among others. For
each job, only one operation can be run at a given moment of time (i.e., jobs
cannot be run “in parallel”), and each machine can be occupied with only one
operation at a time. The timespans $p_{*}$ required for each operation to be
completed are given. The open shop scheduling problem (OSSP) is to find the
schedule with the shortest makespan, i.e., the time interval between the start
of work and the end of the last completed job. For $3$ or more machines and
$3$ or more jobs, OSSP is known to be NP-hard [32].

We have generated a problem set proposed by Taillard [28]222The generator
(Pascal code) is available, e.g., at
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/openshop.txt; upper and
lower bounds for the makespans are provided. For default conversion to CNF, we
have chosen the Crawford-Baker encoding [5].

Specifically, the Crawford-Baker encoding enumerates all operations (a batch
of work for a single job that has to be processed on a certain machine) and
introduces the following Boolean _variables_ :

  * •

$sa_{i,t}$: “op $i$ starts at time $t$ or later”,

  * •

$eb_{i,t}$: “op $i$ ends by time $t$”,

  * •

$pr_{i,j}$: “op $i$ precedes $j$”.

The _constraints_ are formulated as follows:

  * •

$\lor_{t}sa_{i,t}$ for all operations $i$ (each operation begins);

  * •

$\lor_{t}eb_{i,t},~{}t>0$ for all operations $i$ (each operation ends);

  * •

$pr_{i,j}\lor pr_{j,i}$ for all operations $i$, $j$ from the same job that are
run on the same machine (they cannot run at the same time);

  * •

$sa_{i,t}\to sa_{i,t-1}$ (an implication $a\to b$ is the disjunct $\lnot a\lor
b$) for all operations and all times $t$: if $i$ starts at or after $t$, it
starts at or after $t-1$ ($sa$ coherence);

  * •

$eb_{i,t}\to eb_{i,t+1}$ for all possible $i$ and $t$ ($eb$ coherence);

  * •

$sa_{i,t}\to\lnot eb_{i,t+p_{i}-1}$ for all $i$ and $t$, where $p_{i}$ is the
processing time (operation $i$ that started at time $t$ or later cannot end
before $t+p_{i}$);

  * •

$sa_{i,t}\land pr_{i,j}\to sa_{j,t+p_{i}}$ for all $i$, $j$, and times $t$
(operation $j$ cannot start until operation $i$ ends if $i$ precedes $j$).

###  4.3 Graph-Q-SAT implementation modifications

We have made several modifications and improvements to the original _Graph-Q-
SAT_ implementation. First, the Graph-Q-SAT paper [13] makes comparisons
against not the original _Minisat_ solver but against MinisatAgent, an RL
environment agent acting as _Minisat_. The number of steps remains the same,
but the running time required for obtaining the results via MinisatAgent is
significantly increased compared to the original _Minisat_ solver.

To make our approach more feasible for real-life scenarios, we have added a
new action to the environment (trainable in the _CoPilot_ setting); when this
action is triggered, the environment runs the original _Minisat_ without
rebuilding the graph or any calls to _Python_ code whatsoever, which reduces
runtimes very significantly.

In the case of choosing several sequential actions at once, the original
_Graph-Q-SAT_ implementation of both agent and environment resulted in the
most time-consuming stage of each of the agent-environment interaction to
become graph construction. But when we apply several predefined actions one by
one, graph construction for every action is no longer necessary, and actions
can be simply pulled from the action pool; thus, we made the action pool a
part of the environment to avoid unnecessary computations.

##  5 Evaluation results

Table 1: Evaluation on standard SAT benchmarks. | # of | Q-Activity | $\mathbf{SR}(300)$ | $\mathbf{SR}(500)$ | Industrial
---|---|---|---|---|---
| model steps |  | Time | Steps | Time | Steps | Time | Steps
MiniSAT |  |  | 0.069 | 655.1 | 1.292 | 14317.4 | 0.386 | 1033.9
GQSAT | 1 |  | 0.096 | 656.9 | 1.270 | 13770.4 | 0.550 | 1220.2
| 2 |  | 0.121 | 671.0 | 1.377 | 14294.4 | 0.516 | 991.0
| 3 |  | 0.145 | 680.4 | 1.371 | 13918.8 | 0.597 | 1030.8
| 5 |  | 0.180 | 640.7 | 1.429 | 13856.2 | 0.745 | 1070.8
GQSAT-2 | 1 |  | 0.119 | 644.6 | 1.425 | 14804.4 | 0.512 | 981.4
| 2 |  | 0.160 | 645.6 | 1.486 | 14642.7 | 0.628 | 882.9
| 3 |  | 0.203 | 654.7 | 1.693 | 15791.8 | 0.798 | 1074.3
| 5 |  | 0.278 | 661.3 | 1.511 | 13099.7 | 1.097 | 1318.8
CoPilot | 1 |  | 0.121 | 644.6 | 1.422 | 14804.4 | 0.558 | 1177.3
| 2 |  | 0.160 | 645.6 | 1.487 | 14642.7 | 0.639 | 986.9
| 3 |  | 0.201 | 654.7 | 1.708 | 15791.8 | 0.794 | 1026.7
| 5 |  | 0.235 | 654.7 | 1.750 | 15791.8 | 0.904 | 1210.0
CoPilot | 1 | $\checkmark$ | 0.123 | 644.6 | 1.419 | 14804.4 | 0.521 | 1019.3
| 2 | $\checkmark$ | 0.160 | 645.6 | 1.477 | 14642.7 | 0.630 | 926.9
| 3 | $\checkmark$ | 0.204 | 654.7 | 1.704 | 15791.8 | 0.797 | 1102.3
| 5 | $\checkmark$ | 0.231 | 638.8 | 1.516 | 13759.8 | 0.973 | 1217.2
Table 2: Evaluation results for the action pool modification. | # of | Action | $\mathbf{SR}(300)$ | $\mathbf{SR}(500)$ | Industrial
---|---|---|---|---|---
| model | pool |  |  |  |  |  |
| runs | size | Time | Steps | Time | Steps | Time | Steps
MiniSAT |  |  | 0.071 | 759.7 | 0.941 | 11292.7 | 1.881 | 4825.8
GQSAT | 1 |  | 0.088 | 694.1 | 0.944 | 10443.8 | 1.867 | 4828.0
| 2 |  | 0.124 | 611.6 | 0.848 | 9719.1 | 1.931 | 4805.6
| 3 |  | 0.149 | 588.5 | 0.853 | 9507.9 | 2.166 | 4810.9
GQSAT-multistep | 1 | 20 | 0.091 | 611.8 | 0.724 | 8718.5 | 1.842 | 4813.9
| 2 | 20 | 0.107 | 571.4 | 0.721 | 8398.3 | 1.898 | 4532.0
| 3 | 20 | 0.131 | 601.7 | 0.735 | 8409.9 | 1.983 | 4496.9
GQSAT-multistep | 1 | 30 | 0.089 | 587.4 | 0.715 | 8508.3 | 1.771 | 4668.8
| 2 | 30 | 0.110 | 600.6 | 0.732 | 8577.3 | 1.976 | 4734.7
| 3 | 30 | 0.130 | 604.8 | 0.763 | 8681.3 | 3.332 | 7177.7

###  5.1 Results on standard SAT benchmarks

Table 1 shows the main results of our experimental evaluation. In Table 1:

  * •

GQSAT is the original _Graph-Q-SAT_ that uses the RL agent for $1$, $2$, $3$,
or $5$ steps;

  * •

GQSAT-2 denotes our extended _Graph-Q-SAT_ architecture; both versions use our
reimplementation of _Graph-Q-SAT_ as described in Section 4.3;

  * •

MiniSAT denotes the GQSAT implementation where the RL model is not used;

  * •

CoPilot is the GQSAT-2 model with an additional action for releasing control
to MiniSAT, still with a bounded number of steps.

In all cases, MiniSAT is run with restarts; we report average running times
and numbers of iterations. The table shows that restricting the number of RL-
based steps can improve the number of steps while saving time compared to
running it for a longer time; CoPilot can improve the results further.

Another important idea here is that even if we release control to MiniSAT, we
have already spent time on constructing the model and getting its predictions
on first steps. We can reuse the estimated $Q$ functions in MiniSAT’s activity
scores, thus helping it with cold start. For each variable $x_{i}$, we use
$-\frac{1}{\max\left(Q(x_{i}),Q(\lnot x_{i})\right)}$ as the initial value of
its activity score, and then they are modified according to usual MiniSAT
rules. The column “Q-activity” shows this variation of the model, yielding
significant improvements.

Table 3: Evaluation results for the new OSSP representation. | # steps |  $5\times 5$ OSSP |  $6\times 6$ OSSP |  $7\times 7$ OSSP
---|---|---|---|---
MiniSAT |  | 0.886 | 167.4 | 4.483 | 402.7 | 14.444 | 1153.6
GQSAT | 1 | 1.557 | 165.8 | 5.643 | 406.5 | 16.978 | 1152.1
| 2 | 2.174 | 169.7 | 7.201 | 420.5 | 19.028 | 1149.0
| 3 | 2.787 | 164.6 | 8.926 | 431.1 | 22.281 | 1144.2
CoPilot | 1 | 0.890 | 154.7 | 4.423 | 406.8 | 14.423 | 1151.5
with the new | 2 | 1.080 | 170.2 | 4.677 | 403.4 | 14.869 | 1157.0
graph representation | 3 | 1.254 | 171.0 | 5.303 | 387.3 | 16.350 | 1138.9

###  5.2 ML for SAT with an action pool

Table 2 presents the results for our _Graph-Q-SAT_ modification proposed in
Section 3.2, where the agent chooses several actions at once (Fig. 3b). We
show several values of the action pool size (number of “precomputed” actions)
and numbers of model runs (phrased in this way because now a model run
corresponds to several steps); MiniSAT is run without restarts in this table.
Figure 4 shows sample more detailed results for the SR(500) benchmark family
in terms of both number of model runs (Fig. 4a) and running time (Fig. 4b)
shown as a function of the action pool size.

Table 2 shows that the proposed approach with the action pool combines a
significant reduction in the number of steps (much more so than in Table 1)
with virtually no loss in evaluation runtime; we note here that the neural
networks involved in our experiments were run in their default _PyTorch_
implementations without any optimization tricks, and in an industrial
environment ML-assisted runtimes might become significantly faster.

###  5.3 OSSP modification

The OSSP dataset highlights the main weakness of _Graph-Q-SAT_ : even the
smallest nontrivial (i.e., non-polynomial in the worst case) set, _ossp3-3_ ,
has thousands of variables and clauses, which is not a problem for classical
SAT solvers but bloats the _Graph-Q-SAT_ model size and makes each RL agent
step extremely expensive even during inference, let alone training the agent.

Hence, we propose a novel approach to solving OSSP via SAT with GNN-based
heuristics. We still encode OSSP as SAT instances and use MiniSAT with neural
heuristics, but we introduce a more compact graph representation for OSSP.

|
---|---
(a) Number of steps | (b) Running time, s
Figure 4: Experimental results for the SR(500) benchmark.

In our representation, the vertices correspond to operations, and every vertex
is labeled with operation time, first possible timestamp when it can begin,
and last timestamp when it has to finish (deadline). The first number is
constant while the second and third depend on the current solver state; the
labels are initialized with $(p_{i},1,T)$ and then are computed from the
solver state on every step. We add an edge from vertex $i$ to $j$ if operation
$i$ can be done before $j$ and $i$ and $j$ must be done on the same machine or
belong to the same job (add edges between vertices whose operations cannot be
done at the same time). A variable assignment now corresponds to choosing an
edge: an edge from $i$ to $j$ means that operation $i$ must be performed
before op $j$, i.e., we set $pr_{ij}=\mathrm{True}$.

The resulting graph is much smaller than the SAT graph for the Crawford-Baker
encoding: for $j$ jobs and $m$ machines it has $jm$ vertices and $O(jm(j+m))$
edges (note that this size does not depend on $T$) while the SAT formula has
$O(j^{2}m^{2}+jmT)$ variables and $O(j^{2}m^{2}(j+m)+T(jm)^{2})$ clauses. For
instance, a sample $7\times 7$ OSSP task is encoded here with $49$ vertices
and $84$ edges, while the SAT formula graph contains over $54$K variables and
$1.1$M clauses. In fact, _Graph-Q-SAT_ graphs for realistic sized problems
become so large that they do not fit into GPU memory for available hardware.

Table 3 shows that with the new graph representation, CoPilot operates almost
as quickly as the original MiniSAT while reducing the number of iterations and
much faster than _Graph-Q-SAT_ ; note that the difference is significant even
though there were only at most $3$ steps made with the model.

##  6 Conclusion

In this work, we have proposed an approach that limits heavy RL-based
heuristics and finds the best tradeoffs between the number of iterations and
their runtimes for SAT solvers. We have also introduced a novel graph
representation for OSSP-based problems that greatly reduces model size and
speeds up computations. In further work, we plan to extend this exciting
direction, developing similar modifications for other optimization problems
that can be reduced to Boolean satisfiability.

### Acknowledgements

This research has been supported by the Huawei project TC20211214628.

## References

  * [1] S. Amizadeh, S. Matusevych, and M. Weimer.  Learning to solve Circuit-SAT: An unsupervised differentiable approach.  In International Conference on Learning Representations, 2018.
  * [2] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. F. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, Ç. Gülçehre, H. F. Song, A. J. Ballard, J. Gilmer, G. E. Dahl, A. Vaswani, K. R. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. M. Botvinick, O. Vinyals, Y. Li, and R. Pascanu.  Relational inductive biases, deep learning, and graph networks.  CoRR, abs/1806.01261, 2018.
  * [3] A. Biere, K. Fazekas, M. Fleury, and M. Heisinger.  CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling entering the SAT Competition 2020.  In T. Balyo, N. Froleyks, M. Heule, M. Iser, M. Järvisalo, and M. Suda, editors, Proc. of SAT Competition 2020, volume B-2020-1 of Department of Computer Science Report Series B, pages 51–53, 2020.
  * [4] A. Biere and A. Fröhlich.  Evaluating CDCL variable scoring schemes.  In M. Heule and S. Weaver, editors, Theory and Applications of Satisfiability Testing – SAT 2015, pages 405–422, Cham, 2015. Springer International Publishing.
  * [5] J. M. Crawford and A. B. Baker.  Experimental results on the application of satisfiability algorithms to scheduling problems.  In Proceedings of the Twelfth AAAI National Conference on Artificial Intelligence, pages 1092–1097, 1994.
  * [6] D. Devlin and B. O’Sullivan.  Satisfiability as a classification problem.  In Proc. of the 19th Irish Conf. on Artificial Intelligence and Cognitive Science, 2008.
  * [7] N. Eén and N. Sörensson.  An extensible SAT-solver.  In E. Giunchiglia and A. Tacchella, editors, Theory and Applications of Satisfiability Testing, pages 502–518, Berlin, Heidelberg, 2004\. Springer Berlin Heidelberg.
  * [8] J. M. Han.  Enhancing SAT solvers with glue variable predictions.  arXiv:2007.02559, 2020.
  * [9] J. M. Han.  Learning cubing heuristics for SAT from DRAT proofs.  In Conference on Artificial Intelligence and Theorem Proving (AITP), 2020.
  * [10] S. B. Holden.  Machine learning for automated theorem proving: Learning to solve sat and qsat.  Found. Trends Mach. Learn., 14(6):807–989, nov 2021.
  * [11] H. Hoos and T. Stützle.  SATLIB: An online resource for research on SAT, pages 283–292.  04 2000.
  * [12] S. Jaszczur, M. Luszczyk, and H. Michalewski.  Neural heuristics for SAT solving.  CoRR, abs/2005.13406, 2020.
  * [13] V. Kurin, S. Godil, S. Whiteson, and B. Catanzaro.  Can Q-learning with graph networks learn a generalizable branching heuristic for a SAT solver?  Advances in Neural Information Processing Systems, 33:9608–9621, 2020.
  * [14] M. G. Lagoudakis and M. L. Littman.  Learning to select branching rules in the DPLL procedure for satisfiability.  Electronic Notes in Discrete Mathematics, 9:344–359, 2001.
  * [15] J. Liang, V. Ganesh, P. Poupart, and K. Czarnecki.  Exponential recency weighted average branching heuristic for SAT solvers.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 30(1), 2016.
  * [16] J. H. Liang, V. Ganesh, P. Poupart, and K. Czarnecki.  Learning rate based branching heuristic for SAT solvers.  In International Conference on Theory and Applications of Satisfiability Testing, pages 123–140. Springer, 2016.
  * [17] M. Liu, F. Jia, P. Huang, F. Zhang, Y. Sun, S. Cai, F. Ma, and J. Zhang.  Can graph neural networks learn to solve MaxSAT problem?  CoRR, abs/2111.07568, 2021.
  * [18] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. A. Riedmiller.  Playing Atari with deep reinforcement learning.  CoRR, abs/1312.5602, 2013.
  * [19] 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, pages 530–535, 2001.
  * [20] S. Nejati, J. H. Liang, C. Gebotys, K. Czarnecki, and V. Ganesh.  Adaptive restart and CEGAR-based solver for inverting cryptographic hash functions.  In Working Conference on Verified Software: Theories, Tools, and Experiments, pages 120–131. Springer, 2017.
  * [21] E. Ozolins, K. Freivalds, A. Draguns, E. Gaile, R. Zakovskis, and S. Kozlovics.  Goal-aware neural SAT solver.  In 2022 International Joint Conference on Neural Networks (IJCNN). IEEE, 2022.
  * [22] J. M. Robson.  The complexity of Go.  In R. E. A. Mason, editor, Information Processing 83, Proceedings of the IFIP 9th World Computer Congress, Paris, France, September 19-23, 1983, pages 413–417. North-Holland/IFIP, 1983.
  * [23] D. Selsam and N. Bjørner.  Guiding high-performance SAT solvers with unsat-core predictions.  In International Conference on Theory and Applications of Satisfiability Testing, pages 336–353. Springer, 2019.
  * [24] D. Selsam, M. Lamm, B. Bünz, P. Liang, L. de Moura, and D. L. Dill.  Learning a SAT solver from single-bit supervision.  In 7th International Conference on Learning Representations, ICLR 2019, 2019.
  * [25] J. P. M. Silva, I. Lynce, and S. Malik.  Conflict-driven clause learning SAT solvers.  In A. Biere, M. Heule, H. van Maaren, and T. Walsh, editors, Handbook of Satisfiability, volume 185 of Frontiers in Artificial Intelligence and Applications, pages 131–153. IOS Press, 2009.
  * [26] J. A. Storer.  On the complexity of chess.  Journal of Computer and System Sciences, 27(1):77–100, 1983.
  * [27] R. S. Sutton and A. G. Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [28] E. Taillard.  Benchmarks for basic scheduling problems.  European Journal of Operational Research, 64(2):278–285, 1993.
  * [29] H. van Hasselt, A. Guez, and D. Silver.  Deep reinforcement learning with Double Q-learning.  CoRR, abs/1509.06461, 2015.
  * [30] F. Wang and T. Rompf.  From gameplay to symbolic reasoning: Learning SAT solver heuristics in the style of Alpha(Go) Zero.  arXiv:1802.05340, 2018.
  * [31] W. Wang, Y. Hu, M. Tiwari, S. Khurshid, K. L. McMillan, and R. Miikkulainen.  NeuroComb: Improving SAT solving with graph neural networks.  CoRR, abs/2110.14053, 2021.
  * [32] D. P. Williamson, L. A. Hall, J. A. Hoogeveen, C. A. J. Hurkens, J. K. Lenstra, S. V. Sevast’janov, and D. B. Shmoys.  Short shop schedules.  Operations Research, 45(2):288–294, 1997.
  * [33] H. Wu.  Improving SAT-solving with machine learning.  In Proceedings of the 2017 ACM SIGCSE Technical Symposium on Computer Science Education, pages 787–788, 2017.
  * [34] E. Yolcu and B. Póczos.  Learning local search heuristics for Boolean satisfiability.  Advances in Neural Information Processing Systems, 32, 2019.
  * [35] W. Zhang, Z. Sun, Q. Zhu, G. Li, S. Cai, Y. Xiong, and L. Zhang.  NLocalSAT: Boosting local search with solution prediction.  In Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence. International Joint Conferences on Artificial Intelligence Organization, jul 2020.