# Constraining viscous dark energy equation of state in $f(R,L_{m})$ gravity

Lakhan V. Jaybhaye 0000-0003-1497-276X<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Raja Solanki 0000-0001-8849-7688
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Sanjay Mandal 0000-0003-2570-2335
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P.K. Sahoo 0000-0003-2130-8832
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute
of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.

###### Abstract

In this article, we attempt to describe cosmic late-time acceleration of the
universe in the framework of $f(R,L_{m})$ gravity by using an effective
equation of state when the account is taken of bulk viscosity. We presume a
non-linear $f(R,L_{m})$ functional form, specifically,
$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}$, where $\alpha$ is free model
parameter. We obtain the exact solution of our bulk viscous matter dominated
$f(R,L_{m})$ model, and then we utilize the combined $H(z)+Pantheon+Analysis$
data sets to estimate the best fit values of the free parameters of our model.
Then we characterize the behaviour of the matter-energy density, effective
pressure, and the equation of state (EoS) parameter incorporating the viscous
type fluid. The evolution profile of the effective EoS parameter depicts an
acceleration phase of the cosmic expansion whereas the pressure with the
effect of viscosity exhibits negative behaviour that can lead to the
accelerating expansion of the universe. Moreover, the cosmic matter-energy
density shows expected positive behaviour. Further, we investigate the
behaviour of statefinder parameters for the assumed $f(R,L_{m})$ model. We
find that the evolutionary trajectory of the given model lies in the
quintessence region. In addition, we employ the Om diagnostic test that
indicates our model exhibits quintessence behavior. Lastly, we check the
energy condition criteria and find that violation of SEC occurs in the past,
whereas NEC and DEC satisfies the positivity criteria. We find that our
$f(R,L_{m})$ cosmological model with the effect of bulk viscosity provides a
good fit of the recent observational data and can efficiently describe the
cosmic expansion scenario.

Keywords: $f(R,L_{m})$ gravity, bulk viscosity, Equation of state parameter,
Observational datasets, statefinder parameter, Om diagnostic

##  I Introduction

Cosmology has faced a dramatic change when the observational evidence from
type Ia supernovae searches [1, 2] confirmed the accelerating behaviour of the
cosmic expansion. In the last two decades, a plethora of observational results
such as Large Scale Structure [3], Wilkinson Microwave Anisotropy Probe [4],
Cosmic Microwave Background Radiation [5, 6], and the Baryonic Acoustic
Oscillations [7, 8] agrees with the observed cosmic acceleration. The
prominent explanation to describe this accelerating scenario is the presence
of a dark energy component characterized by an equation of state
$\omega=-1.018\pm 0.057$ for a flat universe [9]. Another promising way to
describe the accelerating expansion of the universe by bypassing the
undetected dark energy component is to consider that the more generic action
describes the gravitational field. The cosmological models in which the
Einstein-Hilbert action of general relativity (GR) is modified by introducing
the generic function $f(R)$, where $R$ denotes the Ricci scalar curvature,
first proposed in [10, 11, 12]. The $f(R)$ gravity model is capable to
describe the expansion mechanism without invoking any exotic dark energy
component [13, 14]. Observational signatures of $f(R)$ gravity models along
with the solar system and equivalence principle constraints, presented in the
references [15, 16, Sean, 17]. In the context of solar system tests, the
viable cosmological models of $f(R)$ gravity do exists [18, 19, 20]. Odintsov
et al. have analysed the $H_{0}$ tension and the role of energy conditions in
$f(R)$ gravity models [21, 22]. One can follow the references [23, 24, 25, 26,
27, 28, 29] to see the various implications of cosmological models of $f(R)$
gravity.

A generalization of the curvature based $f(R)$ gravity that incorporates an
explicit coupling of the generic function $f(R)$ with the matter Lagrangian
density $L_{m}$ appeared in [30]. This coupling case was further extended to
the case of arbitrary matter geometry couplings [31]. Harko and Lobo
investigated the curvature-matter couplings in modified gravity from linear
aspects to conformally invariant theories [32]. Models with non-minimal matter
geometry couplings have great astrophysical and cosmological implications.
Harko studied the galactic rotation curves, the matter Lagrangian and the
energy momentum tensor, thermodynamical features, coupling matter and
curvature in Weyl geometry, in the context of non-minimal couplings [33, 34,
35, 36]. Moreover, Bertolami et al. [37] investigated curvature-matter
couplings in modified gravity and Faraoni examined the viability criterion for
modified gravity with an extra force [38]. Further, Harko and Lobo recently
proposed [39] $f(R,L_{m})$ gravity theory that is a generalization matter
curvature coupling theories, where $f(R,L_{m})$ is a generic function that
depends on the Ricci scalar $R$ and the matter Lagrangian $L_{m}$. In this
theory, the covariant divergence of the stress-energy tensor does not
vanishes, an extra force orthogonal to four velocities arises, and the motion
of test particle is non-geodesic. The models of $f(R,L_{m})$ gravity theory
disobey the equivalence principle, and that is restrained by the solar system
experimental tests [40, 41]. Recently, several interesting results on
$f(R,L_{m})$ gravity have been appeared, for instance, see references [42, 43,
44, 45].

In the presented manuscript, we are going to explore the cosmological
$f(R,L_{m})$ model that exhibits viscous type fluid. The introduction of the
coefficient of viscosity in the models of cosmology has a long history. From a
hydrodynamicist’s point of view, there are two viscosity coefficients commonly
appeared in the literature, namely the bulk viscosity coefficient $\zeta$ and
the shear viscosity coefficient $\eta$. By assuming the observationally
supported spatial isotropy of the cosmos, shear viscosity can be omitted.
Whenever a system gets deviated from its thermal equilibrium, then to recover
its thermal equilibrium state an effective pressure is generated. Bulk
viscosity in a cosmological fluid is the manifestation of such an effective
pressure. The idea is to consider the bulk viscosity coefficient $\zeta$ in
the $f(R,L_{m})$ gravity model. We assume that the coefficient of bulk
viscosity $\zeta$ satisfies a scaling law and that reduces the Einstein case
to a form proportional to the Hubble parameter. It has been appeared that this
scaling law is quite useful. One can check the references to review some
interesting viscous fluid cosmological models [46, 47, 48, 49, 50, 51, 52,
53].

The present manuscript is organized in the following manner. In Sec II, we
present the action and basic formulation governing the dynamics in
$f(R,L_{m})$ gravity. In Sec III, we present the Friedmann like equations
corresponding to the flat FLRW universe. In Sec IV, we assume a $f(R,L_{m})$
functional and then we calculate the expression for the Hubble parameter and
the equation of state (EoS) parameter relating the pressure term of bulk
viscous matter with its energy density. In section Sec V, we estimate the
values of the $H_{0}$ and model parameters that obeys with observations, by
incorporating the combined H(z)+Pantheon+Analysis data sets. In addition, we
characterize the behavior of various parameters such as density, effective
pressure, and the EoS parameter. Further in Sec VI, we investigate the $r-s$
parameter trajectory of our $f(R,L_{m})$ model to check the dark energy
behavior recognized by the assumed model. Moreover, in sec VII and sec VIII,
we employ the Om diagnostic test and energy condition criteria. Finally, in
Sec IX, we conclude our findings.

##  II  $f(R,L_{m})$ Gravity Theory

The generic action for $f(R,L_{m})$ gravity read as

$S=\int{f(R,L_{m})\sqrt{-g}d^{4}x}$ (1)

Here $R$ represents the Ricci scalar curvature and $L_{m}$ denotes the matter
Lagrangian .

One can obtained the Ricci scalar $R$ by contracting the Ricci tensor
$R_{\mu\nu}$ as

$R=g^{\mu\nu}R_{\mu\nu}$ (2)

where the Ricci tensor is given by

$R_{\mu\nu}=\partial_{\lambda}\Gamma^{\lambda}_{\mu\nu}-\partial_{\mu}\Gamma^{\lambda}_{\lambda\nu}+\Gamma^{\lambda}_{\mu\nu}\Gamma^{\sigma}_{\sigma\lambda}-\Gamma^{\lambda}_{\nu\sigma}\Gamma^{\sigma}_{\mu\lambda}$
(3)

with $\Gamma^{\alpha}_{\beta\gamma}$ representing the components of Levi-
Civita connection.

Now we obtained the following field equation governing the dynamics of
gravitational interactions, by varying the action (1) with respect to the
metric tensor $g_{\mu\nu}$,

$f_{R}R_{\mu\nu}+(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu})f_{R}-\frac{1}{2}(f-f_{L_{m}}L_{m})g_{\mu\nu}=\frac{1}{2}f_{L_{m}}T_{\mu\nu}$
(4)

Here $f_{R}\equiv\frac{\partial f}{\partial R}$,
$f_{L_{m}}\equiv\frac{\partial f}{\partial L_{m}}$, and $T_{\mu\nu}$
represents the stress-energy tensor for the cosmic fluid, defined by

$T_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}$ (5)

The connection among the energy-momentum scalar $T$, the matter Lagrangian
term $L_{m}$, and the Ricci scalar curvature $R$ acquired by contracting the
field equation (4) as

$Rf_{R}+3\square f_{R}-2(f-f_{L_{m}}L_{m})=\frac{1}{2}f_{L_{m}}T$ (6)

Here $\square
F=\frac{1}{\sqrt{-g}}\partial_{\alpha}(\sqrt{-g}g^{\alpha\beta}\partial_{\beta}F)$
for any scalar function $F$ .

In addition, one can obtain the following relation by employing the covariant
derivative in equation (4)

$\nabla^{\mu}T_{\mu\nu}=2\nabla^{\mu}ln(f_{L_{m}})\frac{\partial
L_{m}}{\partial g^{\mu\nu}}$ (7)

##  III Motion equations in $f(R,L_{m})$ gravity

In order to probe the cosmological implications, we consider the following
homogeneous and spatially isotropic FLRW metric [54]

$ds^{2}=-dt^{2}+a^{2}(t)[dx^{2}+dy^{2}+dz^{2}]$ (8)

where, $a(t)$ is the cosmic scale factor. The Ricci scalar obtained for the
metric (8) is

$R=6\frac{\ddot{a}}{a}+6\bigl{(}\frac{\dot{a}}{a}\bigr{)}^{2}=6(\dot{H}+2H^{2})$
(9)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter.

The energy-momentum tensor comprises of energy density $\rho$ and the pressure
$\bar{p}$ of the cosmic fluid with viscosity effect is given by,

$\mathcal{T}_{\mu\nu}=(\rho+\bar{p})u_{\mu}u_{\nu}+\bar{p}g_{\mu\nu}$ (10)

where $\bar{p}=p-3\zeta H$ and $u^{\mu}=(1,0,0,0)$ are components of the four
velocities. Here $p$ is the usual pressure and $\zeta>0$ is the coefficient of
bulk viscosity.

The connection between matter-energy density and the usual pressure is given
as [55]

$p=(\gamma-1)\rho$ (11)

where $\gamma$ is a constant with $0\leq\gamma\leq 2$. Hence the effective
equation of state characterizing the bulk viscous cosmic fluid reads as [56,
57, 58]

$\bar{p}=(\gamma-1)\rho-3\zeta H$ (12)

Under the constraint of homogeneity and spatial isotropy, the cosmic fluid
incorporating viscosity possesses dissipative phenomenon. Considering
viscosity in a cosmic fluid can minimize the ideal characteristics of a fluid,
and participates to the total pressure negatively. This can be checked in
References [59, 60, 61].

The Friedmann equations that characterizes the bulk viscous matter dominated
universe in $f(R,L_{m})$ gravity reads as [62]

$3H^{2}f_{R}+\frac{1}{2}\left(f-f_{R}R-f_{L_{m}}L_{m}\right)+3H\dot{f_{R}}=\frac{1}{2}f_{L_{m}}\rho$
(13)

and

$\dot{H}f_{R}+3H^{2}f_{R}-\ddot{f_{R}}-3H\dot{f_{R}}+\frac{1}{2}\left(f_{L_{m}}L_{m}-f\right)=\frac{1}{2}f_{L_{m}}\bar{p}$
(14)

##  IV Cosmological $f(R,L_{m})$ Model

We choose the following $f(R,L_{m})$ function in order to explore the dynamics
of the universe possesses viscosity [62, 63],

$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}$ (15)

Here $\alpha$ is free model parameter. The model under consideration is more
general in nature and it is motivated by the functional form
$f(R,L_{m})=f_{1}(R)+f_{2}(R)G(L_{m})$ that represents arbitrary matter-
geometry coupling [63].

Then for this specific functional form with $L_{m}=\rho$ [64], the Friedmann
equations (13) and (14) characterizing the universe dominated with bulk
viscous matter becomes

$3H^{2}=(2\alpha-1)\rho^{\alpha}$ (16)

and

$2\dot{H}+3H^{2}=\left\\{(\alpha-1)\rho-\alpha\bar{p}\right\\}\rho^{\alpha-1}$
(17)

Now by using equation (7), we obtained the following matter conservation
equation for our bulk viscous cosmological $f(R,L_{m})$ model

$(2\alpha-1)\dot{\rho}+3\gamma H\rho=0$ (18)

From equations (16) and (17), one can have

$\dot{H}+\frac{3\alpha\gamma}{2(2\alpha-1)}H^{2}=\frac{3}{2}\left(\frac{3}{2\alpha-1}\right)^{\frac{\alpha-1}{\alpha}}\alpha\zeta
H^{\frac{3\alpha-2}{\alpha}}$ (19)

We substitute $\frac{1}{H}\frac{d}{dt}=\frac{d}{dln(a)}$ so that equation (19)
becomes

$\frac{dH}{dln(a)}+\frac{3\alpha\gamma}{2(2\alpha-1)}H=\frac{3}{2}\left(\frac{3}{2\alpha-1}\right)^{\frac{\alpha-1}{\alpha}}\alpha\zeta
H^{\frac{2(\alpha-1)}{\alpha}}$ (20)

On integrating the equation (20) we obtained the expression for Hubble
parameter as follows

$H(z)=\big{\\{}H_{0}^{\frac{2-\alpha}{\alpha}}(1+z)^{\frac{3\gamma(2-\alpha)}{2(2\alpha-1)}}+\frac{3\zeta}{\gamma}\left(\frac{2\alpha-1}{3}\right)^{\frac{1}{\alpha}}[1-(1+z)^{\frac{3\gamma(2-\alpha)}{2(2\alpha-1)}}]\big{\\}}^{\frac{\alpha}{2-\alpha}}$
(21)

where $H(0)=H_{0}$ represents the present value of the Hubble parameter. In
particular, for the case $\alpha=1$ with $\gamma=1$ and $\zeta=0$, the
solution reduces to $H(z)=H_{0}(1+z)^{\frac{3}{2}}$ , the usual ordinary
matter dominated universe.

The effective equation of state parameter for our bulk viscous cosmological
model is given by

$\omega_{eff}=\frac{p_{eff}}{\rho}=\gamma-1-\frac{3\zeta H}{\rho}$ (22)

By using equations (16) and (21), one can acquired

$\omega_{eff}=\gamma-1-3\zeta\left(\frac{2\alpha-1}{3}\right)^{\frac{1}{\alpha}}\big{\\{}H_{0}^{\frac{2-\alpha}{\alpha}}(1+z)^{\frac{3\gamma(2-\alpha)}{2(2\alpha-1)}}+\frac{3\zeta}{\gamma}\left(\frac{2\alpha-1}{3}\right)^{\frac{1}{\alpha}}[1-(1+z)^{\frac{3\gamma(2-\alpha)}{2(2\alpha-1)}}]\big{\\}}^{-1}$
(23)

##  V Data, Methodology, and Physical Interpretation

In this section, we estimate the parameter values of our model that is
appropriate to describe the various cosmic epochs, by invoking the $H(z)$ and
Pantheon+Analysis data sets. To calculate the suitable values of $H_{0}$ and
model parameters $\alpha$, $\gamma$, and $\zeta$, we incorporate $31$ points
of $H(z)$ data sets and $1701$ points from the Pantheon+Analysis samples. To
estimate the mean values of the parameters of our viscosity model, we apply
the Bayesian technique and likelihood function along with the Markov Chain
Monte Carlo (MCMC) method in emcee python library [65].

###  V.0.1 H(z) datasets

It is well known that the Hubble parameter can directly investigate cosmic
expansion. In terms of redshift the Hubble parameter can be acquired as
$H(z)=-\frac{1}{1+z}\frac{dz}{dt}$. Since $dz$ is derived from the
spectroscopic surveys therefore one can obtain the model-independent $H(z)$
value by measuring the $dt$. In this manuscript, we incorporate $31$ data
points of $H(z)$ measurements in the redshift range $0.07\leq z\leq 2.41$
[66]. One can check the reference [67] for the complete list of $31$ data
points. We define the chi-square function to find out the mean values of the
bulk viscous model parameters $\alpha$, $\gamma$, $\zeta$, and $H_{0}$ as
follows,

$\chi_{H}^{2}(H_{0},\alpha,\gamma,\zeta)=\sum\limits_{k=1}^{31}\frac{[H_{th}(z_{k},H_{0},\alpha,\gamma,\zeta)-H_{obs}(z_{k})]^{2}}{\sigma_{H(z_{k})}^{2}}.$
(24)

Here, the theoretical value of the $H(z)$ acquired by our cosmological model
is represented by $H_{th}$ whereas $H_{obs}$ denotes its observed value and
$\sigma_{H(z_{k})}$ is the standard error.

###  V.0.2 Pantheon datasets

Earlier, the observational results on type Ia supernovae confirmed that our
universe is going through a phase of accelerated expansion. In the past two
decades, observations on supernovae samples have been extensively increased.
In 2018, $1048$ samples of type Ia supernovae covering the redshift range
$0.01<z<2.3$ has been released which is known as Pantheon supernovae samples
[68]. The PanSTARSS1 Medium, Deep Survey, SDSS, HST surveys, SNLS, and
numerous low redshift surveys contribute to it. Recently, Pantheon+ Analysis
sample incorporating 1701 light curves of $1550$ supernovae in the range of
redshift $[0.001,2.26]$ has been released [69]. The luminosity distance is
taken to be [planck_collaboration/2020],

$\displaystyle D_{L}(z)$ $\displaystyle=$
$\displaystyle\frac{c(1+z)}{H_{0}}S_{k}\left(H_{0}\int_{0}^{z}\frac{1}{H(z^{\prime})}dz^{\prime}\right),$
$\displaystyle\text{where }S_{k}(x)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{c}\sinh(x\sqrt{\Omega_{k}})/\Omega_{k}\text{,
}\Omega_{k}>0\\\ x\text{, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
}\Omega_{k}=0\\\ \sin
x\sqrt{\left|\Omega_{k}\right|})/\left|\Omega_{k}\right|\text{,
}\Omega_{k}<0\end{array}\right.$

For a spatially flat universe, we have

$D_{L}(z)=(1+z)\int_{0}^{z}\frac{cdz^{\prime}}{H(z^{\prime})},$ (26)

where $c$ is the speed of light.

We have calculated the $\chi^{2}$ function for the Pantheon supernovae samples
by correlating the theoretical distance modulus

$\mu(z)=5log_{10}D_{L}(z)+\mu_{0},$ (27)

with

$\mu_{0}=5log(1/H_{0}Mpc)+25,$ (28)

such that

$\chi^{2}_{SN}(p_{1},....)=\sum_{i,j=1}^{1701}\bigtriangledown\mu_{i}\left(C^{-1}_{SN}\right)_{ij}\bigtriangledown\mu_{j},$
(29)

Here $p_{j}$ represents free model parameters and $C_{SN}$ is the covariance
matrix [69], and

$\displaystyle\quad\bigtriangledown\mu_{i}=\mu^{th}(z_{i},p_{1},...)-\mu_{i}^{obs}.$

where $\mu_{th}$ represents value of the distance modulus predicted by our
model while $\mu_{obs}$ its observed value.

Now the $\chi^{2}$ function for the H(z)+Pantheon+Analysis data sets is taken
to be

$\chi^{2}_{total}=\chi^{2}_{H}+\chi^{2}_{SN}$ (30)

We present the $1-\sigma$ and $2-\sigma$ likelihood contours for the model
parameters $\alpha$, $\gamma$, $\zeta$, and $H_{0}$ using combined
H(z)+Pantheon+Analysis data sets below.

Figure 1: The $1-\sigma$ and $2-\sigma$ contours for the model parameters
$\alpha$, $\gamma$, $\zeta$, and $H_{0}$ using combined H(z)+Pantheon+Analysis
data sets

The obtained best fit values are $\alpha=1.310^{+0.037}_{-0.032}$,
$\gamma=1.29\pm 0.20$, $\zeta=5.02\pm 0.26$, and $H_{0}=72.09\pm 0.19$.

Now we are going to present the cosmological implications of obtained
observational constraints. We analyze the behaviour of energy density,
pressure component incorporating viscosity, and the effective EoS parameter
for the obtained mean values of $H_{0}$ and model parameters $\alpha$,
$\gamma$, and $\zeta$ constrained by the H(z)+Pantheon+Analysis data sets.

Figure 2: The reconstruction of the energy density as a function of the
redshift for our model is presented for 7500 samples which are reproduced by
re-sampling the chains through emcee. We plot all the obtained curves,
alongside the curve corresponding to the best fit of the parameters (red
curve). Figure 3: The reconstruction of the effective pressure as a function
of the redshift for our model is presented for 7500 samples which are
reproduced by re-sampling the chains through emcee. We plot all the obtained
curves, alongside the curve corresponding to the best fit of the parameters
(red curve). Figure 4: The reconstruction of the effective EoS parameter as a
function of the redshift for our model is presented for 7500 samples which are
reproduced by re-sampling the chains through emcee. We plot all the obtained
curves, alongside the curve corresponding to the best fit of the parameters
(red curve).

We reconstructed the matter-energy density, the effective pressure and EoS
parameter as a function of the redshift, presented in Figs. 2, 3, and 4, for
7500 samples that are reproduced by re-sampling the chains through emcee. From
Fig. 2 it is evident that the cosmic matter-energy density shows expected
positive behaviour and it vanishes with the expansion of the universe in the
far future. The effective pressure component presented in Fig. 3 exhibits
negative behaviour that can lead to the accelerating expansion of the
universe. Further, the present value of the effective EoS parameter is
obtained to be $\omega_{0}\approx-0.71$. Thus, the behaviour of the effective
EoS parameter in Fig. 4 confirmed the accelerating nature of the expansion
phase of the universe.

##  VI Statefinder Diagnostic

It is well accepted that the responsible behind cosmic expansion is the dark
energy. In the last few decades, the investigation of the origin and
fundamental behavior of dark energy is increased. Consequently, plenty of dark
energy models started appearing, and therefore the either quantitative or
qualitative distinction between these models of dark energy becomes necessary.
In this direction, Sahni et al. [70] proposed a statefinder diagnostic method
that can identify amongst various dark energy models with the help of a pair
of geometrical parameters called statefinder parameters $(r,s)$. It is defined
as

$r=\frac{\dddot{a}}{aH^{3}}$ (31)

and

$s=\frac{(r-1)}{3(q-\frac{1}{2})}$ (32)

We evaluate the statefinder parameters $(r,s)$ for our cosmological
$f(R,L_{m})$ model. The evolutionary trajectory of the assumed model with the
agreement of obtained observational constraints is presented in Fig. 5. The
deviation of the evolutionary trajectory of the given model from the
$\Lambda$CDM one gives the required discrimination. The value $r=1,s=0$
represents the $\Lambda$CDM model, $r>1,s<0$ represents the Chaplygin gas
model, and $r<1,s>0$ represents the quintessence model. The present value of
statefinder parameters for our model is nearly $(r,s)=(0.43,0.33)$. From Fig.
5 it is evident that the dark component due to modified geometry with the
effect of bulk viscosity has quintessence type behavior.

Figure 5: Profile of the evolution trajectory of given model in the $r-s$
plane with the agreement of obtained observational constraints.

##  VII Om Diagnostics

The Om diagnostic is another recently proposed method that can effectively
distinguish different dark energy models [71]. It is much simpler as compared
to that of statefinder analysis since it offers the formulation incorporating
only the Hubble parameter. For spatially flat constraint, it is given by,

$Om(z)=\frac{\big{(}\frac{H(z)}{H_{0}}\big{)}^{2}-1}{(1+z)^{3}-1}$ (33)

The negative slope of $Om(z)$ represents quintessence behaviour whereas
positive slope represents phantom behaviour . The constant nature of $Om(z)$
corresponds to the $\Lambda$CDM type behaviour of the given model.

Figure 6: Profile of Om diagnostic parameter with the agreement of obtained
observational constraints.

Fig.6 indicate that the Om diagnostic parameter shows negative slope in the
entire domain. Thus, Om diagnostic test indicate that our bulk viscous matter
dominated $f(R,L_{m})$ model follows quintessence scenario.

##  VIII Energy Conditions

Now, we are going to test the viability of the acquired solution corresponding
to the assumed $f(R,L_{m})$ model by invoking energy conditions criterion. The
energy conditions are criteria imposed to the energy-momentum tensor, in order
to fulfill the positivity condition of energy. These criteria are offered from
the excellent work of Raychaudhuri that is known as Raychaudhuri’s equation
and are written as [72]

  * •

Null energy condition (NEC) : $\rho_{eff}+p_{eff}\geq 0$;

  * •

Weak energy condition (WEC) : $\rho_{eff}\geq 0$ and $\rho_{eff}+p_{eff}\geq
0$;

  * •

Dominant energy condition (DEC) : $\rho_{eff}\pm p_{eff}\geq 0$;

  * •

Strong energy condition (SEC) : $\rho_{eff}+3p_{eff}\geq 0$,

with $\rho_{eff}$ is the effective energy density.

Figure 7: Behavior of the NEC vs redshift. Figure 8: Behavior of the DEC vs
redshift. Figure 9: Behavior of the SEC vs redshift.

From Figs.7 and 8, we observed that the NEC and DEC satisfies the positivity
criteria in the entire domain of redshift range, corresponding to the
estimated values of parameters from observational data sets. As WEC comprises
of energy density and NEC, it is also satisfied. Finally Fig.9 shows that the
violation of SEC occurs in the recent past, and hence this violations favors
the cosmic acceleration.

##  IX Conclusion

Hydrodynamically, the inclusion of the coefficient of viscosity in the cosmic
matter content is quite natural, as an ideal characteristics of a fluid is
after all an abstraction. In the presented article, we have analyzed the
significance of bulk viscosity to drive the cosmic late time acceleration
under the $f(R,L_{m})$ background. $f(R,L_{m})$ gravity theory is a
generalization of matter curvature coupling theories [39]. Harko and Lobo
investigated the curvature-matter couplings in modified gravity from linear
aspects to conformally invariant theories [32]. Models with non-minimal matter
geometry couplings have great astrophysical and cosmological implications [35,
36, 37, 38]. For our analysis, we considered a $f(R,L_{m})$ function,
particularly, $f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}$, where $\alpha$ is free
model parameter. Then we have assumed the effective equation of state in
equation 12, which is the Einstein case value with proportionality constant
$\zeta$ used in the Einstein theory [46] frequently used in the literature. We
found the exact solution of our bulk viscous matter dominated $f(R,L_{m})$
model, and then we used the combined $H(z)+Pantheon+Analysis$ observational
data sets to constrain the present value of the Hubble parameter $H_{0}$ and
the model parameters. The obtained best fit values are
$\alpha=1.310^{+0.037}_{-0.032}$, $\gamma=1.29\pm 0.20$, $\zeta=5.02\pm 0.26$,
and $H_{0}=72.09\pm 0.19$. In addition, we have characterized the behaviour of
matter-energy density, pressure component incorporating viscosity, and the
effective EoS parameter as a function of the redshift, presented in Figs. 2,
3, and 4, for 7500 samples that are reproduced by re-sampling the chains
through emcee. From Fig. 2 it is evident that the cosmic matter-energy density
shows expected positive behaviour and the effective pressure component
presented in Fig. 3 exhibits negative behaviour that can lead to the
accelerating expansion of the universe. Moreover, the present value of the
effective EoS parameter is obtained to be $\omega_{0}\approx-0.71$. Thus, the
trajectory of the EoS parameter in Fig. 4 confirmed the accelerating nature of
the expansion phase of the universe. Then we evaluated the $(r,s)$ parameters
for our assumed $f(R,L_{m})$ model. The present value of statefinder
parameters for our model is nearly $(r,s)=(0.43,0.33)$. From Fig. 5 we noticed
that the evolutionary trajectory of our $f(R,L_{m})$ model lies in the
quintessence region. Further, the Om diagnostic presented in Fig 6, indicates
that our assumed $f(R,L_{m})$ model favors the quintessence type dark energy.
Finally, the energy conditions presented in Fig 7, 8, and 9, exhibits
positivity criteria in the entire domain of redshift range corresponding to
the case of NEC and DEC, whereas it shows the violation in case of SEC. This
violation of SEC, occurs in the recent past, favors the observed acceleration.
We conclude that our cosmological $f(R,L_{m})$ model with the fluid
incorporating the bulk viscosity effects, can efficiently interpret the late
time cosmic phenomenon of the universe with observational compatibility.

## Data Availability Statement

There are no new data associated with this article.

## Acknowledgments

L.V.J. acknowledges University Grant Commission (UGC), Govt. of India, New
Delhi, for awarding JRF (NTA Ref. No.: 191620024300). R.S. acknowledges UGC,
New Delhi, India for providing Junior Research Fellowship (UGC-Ref. No.:
191620096030). SM acknowledges Department of Science and Technology (DST),
Govt. of India, New Delhi, for awarding Senior Research Fellowship (File No.
DST/INSPIRE Fellowship/2018/IF18D676). We are very much grateful to the
honorable referees and to the editor for the illuminating suggestions that
have significantly improved our work in terms of research quality, and
presentation.

## References

  * [1] Riess, A.G.; Filippenko, A.V.; Challis, P.; Clocchiatti, A.; Diercks, A.; Garnavich, P.M.; Gilliland, R.L.; Hogan, C.J.; Jha, S.; Kirshner, R.P.; et al. Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant. Astron. J. 1998, _116_ , 1009.
  * [2] Perlmutter, S.; Aldering, G.; Goldhaber, G.; Knop, R.A.; Nugent, P.; Castro, P.G.; Deustua, S.; Fabbro, S.; Goobar, A.; Groom, D.E.; et al. Measurements of $\Omega$ and $\Lambda$ from 42 High-Redshift Supernovae. Astrophys. J. 1999, _517_ , 565.
  * [3] Scoville, N.; Aussel, H.; Brusa, M.; Capak, P.; Carollo, C. M.; Elvis, M.; Giavalisco, M.; Guzzo, L.; Hasinger, G.; Impey, C. The Cosmic Evolution Survey (COSMOS): Overview. Astrophys. J., Suppl. Ser. 2007, _172_ , 1.
  * [4] Spergel, D.N.; Verde, L.; Peiris, H.V.; Komatsu, E.; Nolta, M. R.; Bennett, C.L.; Halpern, M.; Hinshaw, G.; Jarosik, N.; Kogut, A. First-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Determination of Cosmological Parameters. Astrophys. J. Suppl. 2003, _148_ , 175.
  * [5] Bennett, C.L.; Halpern, M.; Hinshaw, G.; Jarosik, N.; Kogut, A.; Limon, M.; Meyer, S.S.; Page, L.; Spergel, D.N.; Tucker, G.S. First-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Preliminary Maps and Basic Results Astrophys. J. Suppl. 2003, _148_ , 119–134.
  * [6] Caldwell, R.R.; Doran, M. Cosmic microwave background and supernova constraints on quintessence: Concordance regions and target models. Phys. Rev. D 2004, _69_ , 103517.
  * [7] Eisenstein, D.J.; Zehavi, I.; Hogg, D.W.; Scoccimarro, R.; Blanton, M.R.; Nichol, R.C.; Scranton, R.; Seo, H.J.; Tegmark, M.; Zheng, Z. Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies. Astrophys. J. 2005, _633_ , 560.
  * [8] Percival, W.J.; Reid, B.A.; Eisenstein, D.J.; Bahcall, N.A.; Budavari, T.; Frieman, J.A.; Fukugita, M.; Gunn, J.E.; Ivezic, Z.; Knapp, G.R.; et al. Baryon acoustic oscillations in the Sloan Digital Sky Survey Data Release 7 galaxy sample. Mon. Not. R. Astron. Soc. 2010, _401_ , 2148.
  * [9] Aghanim, N.; Akrami, Y.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Basak, S.; et al. Planck Collaboration. Astron. Astrophys. 2020, _641_ , A6.
  * [10] Buchdahl, H.A. Non-Linear Lagrangians and Cosmological Theory. Mon. Not. R. Astron. Soc. 1970, _150_ , 1.
  * [11] Kerner, R. Cosmology without singularity and nonlinear gravitational Lagrangians. Gen. Relativ. Gravit. 1982, _14_ , 453.
  * [12] Kleinert, H.; Schmidt, H.J. Cosmology with Curvature-Saturated Gravitational Lagrangian $R$. Gen. Relativ. Gravit. 2002, _34_ , 1295.
  * [13] Carroll, S.M.; Duvvuri, V.; Trodden, M.; Turner M.S. Is cosmic speed-up due to new gravitational physics ?. Phys. Rev. D 2004, _70_ , 043528.
  * [14] Capozziello, S.; Nojiri, S.; Odintsov, S.D.; Troisi, A. Cosmological viability of $f(R)$-gravity as an ideal fluid and its compatibility with a matter dominated phase. Phys. Lett. B 2006, _639_ , 135.
  * [15] Tsujikawa, S. Observational signatures of $f(R)$ dark energy models that satisfy cosmological and local gravity constraints. Phys. Rev. D 2008, _77_ , 023507.
  * [16] Capozziello, S.; Tsujikawa, S. Solar system and equivalence principle constraints on gravity by the chameleon approach. Phys. Rev. D 2008, _77_ , 107501.
  * [17] Starobinsky, A.A. Disappearing cosmological constant in $f(R)$ gravity. JETP Lett. 2007, _86_ , 157–163.
  * [18] Nojiri, S.; Odintsov, S.D. Modified gravity with negative and positive powers of curvature: Unification of inflation and cosmic acceleration. Phys. Rev. D 2003, _68_ , 123512.
  * [19] Faraoni, V. Solar system experiments do not yet veto modified gravity models. Phys. Rev. D 2006, _74_ , 023529.
  * [20] Amendola, L.; Tsujikawa, S. Phantom crossing, equation-of-state singularities, and local gravity constraints in $f(R)$ models. Phys. Lett. B 2008, _660_ , 125.
  * [21] Odintsov, S.D.; Gomez, D.S.C.; Sharov, G.S. Analyzing the $H_{0}$ tension in $f(R)$ gravity models. Nucl. Phys. B 2021, _966_ , 115377.
  * [22] Capozziello, S.; Nojiri, S.; Odintsov, S.D. The role of energy conditions in $f(R)$ cosmology. Phys. Lett. B 2018, _781_ , 99–106.
  * [23] Nojiri, S.; Odintsov, S. Unified cosmic history in modified gravity: from $f(R)$ theory to Lorentz non-invariant models.Phys. Lett.B 2007, _657_ , 238.
  * [24] Nojiri, S.; Odintsov, S.D. Modified $f(R)$ gravity unifying $R^{m}$ inflation with the $\Lambda CDM$ epoch. Phys. Rev. D 2008, _77_ , 026007.
  * [25] Santos, J.; Alcaniz,J.S.; Reboucas, M.J.; Carvalho, F. C. Energy conditions in $f(R)$ gravity. Phys. Rev. D 2007, _76_ , 083513.
  * [26] Nojiri, S.; Odintsov, S.; Oikonomou, V.K. Unifying inflation with early and late-time dark energy in $f(R)$ gravity. Phys. Dark Univ. 2020, _29_ , 100602.
  * [27] Odintsov, S.D.; Oikonomou, V.K. Inflationary attractors in $f(R)$ gravity. Phys. Lett. B 2020, _807_ , 135576.
  * [28] Odintsov, S.D.; Oikonomou, V.K., Dark energy oscillations in mimetic $f(R)$ gravity. Phys. Rev. D 2016, _94_ , 044012.
  * [29] Paliathanasis, A. Similarity solutions for the Wheeler-DeWitt equation in $f(R)$-cosmology. Eur. Phys. J. C 2019, _79_ , 1031.
  * [30] Bertolami, O.; Boehmer,C.G.; Harko, T.; Lobo, F.S.N. Extra force in $f(R)$ modified theories of gravity. Phys. Rev. D 2007, _75_ , 104016.
  * [31] Harko, T. Modified gravity with arbitrary coupling between matter and geometry. Phys. Lett. B 2008, _669_ , 376.
  * [32] Lobo, F.S.N.; Harko, T. Palatini formulation of the conformally invariant $f(R,L_{m})$ gravity theory. Int. J. Mod. Phys. D 2022, _31_ , 2240010.
  * [33] Harko, T. Galactic rotation curves in modified gravity with nonminimal coupling between matter and geometry. Phys. Rev. D 2010, _81_ , 084050.
  * [34] Harko, T. The matter Lagrangian and the energy-momentum tensor in modified gravity with nonminimal coupling between matter and geometry. Phys. Rev. D 2010, _81_ , 044021.
  * [35] Harko, T. Thermodynamic interpretation of the generalized gravity models with geometry-matter coupling. Phys. Rev. D 2014, _90_ , 044067.
  * [36] Harko, T.; Shahidi, S. Coupling matter and curvature in Weyl geometry: conformally invariant $f(R,L_{m})$ gravity. Eur. Phys. J. C 2022, _82_ , 219.
  * [37] Bertolami, O.; Harko, T.; Lobo, F.S.N.; Paramos, J. Non-minimal curvature-matter couplings in modified gravity. arXiv, 2008, arXiv:0811.2876v1.
  * [38] Faraoni, V. Viability criterion for modified gravity with an extra force. Phys. Rev. D 2007, _76_ , 127501.
  * [39] Harko, T.; Lobo, F.S.N., $f(R,L_{m})$ gravity. Eur. Phys. J. C 2010, _70_ , 373–379.
  * [40] Faraoni, V. Cosmology in Scalar-Tensor Gravity; Kluwer Academic: Dordrecht, The Netherlands, 2004.
  * [41] Bertolami, O.; PÃ¡ramos, J.; Turyshev, S. General Theory of Relativity: Will it survive the next decade?. arXiv, 2006, arXiv:gr-qc/0602016.
  * [42] Goncalves, B.S.; Moraes, P.H.R.S. Cosmology from non-minimal geometry-matter coupling. arXiv, 2021, arXiv:2101.05918.
  * [43] Labato, R.V.; Carvalho, G. A.; Kelkar, N. G.; Nowakowski, M. Massive white dwarfs in $f(R,L_{m})$ gravity. Eur. Phys. J. C 2022, _82_ , 540.
  * [44] Labato, R.V.; Carvalho, G.A.; Bertulani, C.A., Neutron stars in $f(R,L_{m})$ gravity with realistic equations of state: joint-constrains with GW170817, massive pulsars, and the PSR J0030+0451 mass-radius from NICER data. Eur. Phys. J. C 2021, _81_ , 1013.
  * [45] Jaybhaye, L.V.; S. mandal; Sahoo, P.K., Constraints on Energy Conditions in $f(R,L_{m})$ Gravity. Int. J. Geom. Methods Mod. 2022, _19_ , 2250050.
  * [46] Brevik, I. Viscosity in Modified Gravity. Entropy 2012, _2012_ , 2302–2310.
  * [47] Brevik, I.; Gron, O., Universe Models with Negative Bulk Viscosity. Astrophys. Space Sci. 2013, _347_ , 399.
  * [48] Brevik, I.; Gron, O.; Haro, J.; Odintsov, S.D.; Saridakis, E.N. Viscous cosmology for early- and late-time universe. Int. J. Mod. Phys. D 2017, _26_ , 1730024.
  * [49] Brevik, I.; Makarenko, A.N.; Timoshkin, A.V. Viscous accelerating universe with nonlinear and logarithmic equation of state fluid. Int. J. Geom. Methods Mod. 2019, _16_ , 1950150
  * [50] Brevik, I.; Normann, B.D. Remarks on cosmological bulk viscosity in different epochs. Symmetry 2020, _2020_ , 1085.
  * [51] Mohan, N.D.J.; Sasidharan, A.; Mathew, T.K. Bulk viscous matter and recent acceleration of the universe based on causal viscous theory. Eur. Phys. J. C 2017, _77_ , 849.
  * [52] Astashenok, A.V.; Odintsov, S.D.; Tepliakov, A.S. The unified history of the viscous accelerating universe and phase transitions. Nucl. Phys. B 2022, _974_ , 115646.
  * [53] Sasidharan, A.; Mathew, T.K. Bulk viscous matter and recent acceleration of the universe. Eur. Phys. J. C 2015, _75_ , 348.
  * [54] Ryden, B. Introduction to Cosmology; Addison Wesley: San Francisco, CA, USA, 2003.
  * [55] Ren, J.; Meng, X.H. Cosmological model with viscosity media (dark fluid) described by an effective equation of state. Phys. Lett. B 2006, _633_ , 1–8.
  * [56] Brevik, I.; Gorbunova, O. Dark Energy and Viscous Cosmology. Gen. Rel. Grav. 2005, _37_ , 2039.
  * [57] Gron, O. Viscous inflationary universe models. Astrophys. Sci. 1990, _173_ , 191-225.
  * [58] Eckart, C. The Thermodynamics of Irreversible Processes. III. Relativistic Theory of the Simple Fluid. Phys. Rev. 1940, _58_ , 919.
  * [59] Odintsov, S.D.; Gomez, D.S.C.; Sharov, G.S. Testing the equation of state for viscous dark energy. Phys. Rev. D 2020, _101_ , 044010.
  * [60] Fabris, J.C.; Goncalves, S.V.B.; Ribeiro, R.S. Bulk viscosity driving the acceleration of the Universe. Gen. Rel. Grav. 2006, _3_ , 495.
  * [61] Meng, X.H.; Dou, X. Singularity and entropy of the viscosity dark energy model. Comm. Theor. Phys. 2009, _52_ , 377.
  * [62] Jaybhaye, L.V.; Solanki,R.; Mandal, S.; Sahoo, P.K. Cosmology in $f(R,L_{m})$ gravity. Phys. Lett. B 2022, _831_ , 137148.
  * [63] Harko, T.; Lobo, F.S.N. Generalized Curvature-Matter Couplings in Modified Gravity. Galaxies 2014, _2014_ , 410-465.
  * [64] Harko, T.; Lobo, F.S.N.; Mimoso, J.P.; Pavon, D. Gravitational induced particle production through a nonminimal curvature-matter coupling. Eur. Phys. J. C 2015, _75_ , 386.
  * [65] Mackey, D.F.; HOGG, D.W.; LANG, D.; GOODMAN, J. emcee: The MCMC Hammer. Publ. Astron. Soc. Pac. 2013, _125_ , 306.
  * [66] Sharov, G.S.; Vasiliev, V.O. How predictions of cosmological models depend on Hubble parameter data sets. Math. Model. Geom. 2018, _6_ , 1.
  * [67] Solanki, R.; Pacif, S.K.J.; Parida, A.; Sahoo, P.K. Cosmic acceleration with bulk viscosity in modified $f(Q)$ gravity Phys. Dark Univ. 2021, _32_ , 100820.
  * [68] Scolnic, D.M.; Jones, D.O.; Rest, A.; Pan, Y.C.; Chornock, R.; Foley, R.J.; Huber, M.E.; Kessler, R.; Narayan, G.; Riess, A.G.; et al. The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from PanSTARRS1 and Cosmological Constraints from the Combined Pantheon Sample. Astrophys. J. 2018, _859_ , 101.
  * [69] Brout, D.; Scolnic, D.; Popovic, B.; Riess, A.G.; Carr, A.; Zuntz, J.; Kessler, R.; Davis, T.M.; Hinton, S.; Jones, D.; et al. The Pantheon+ Analysis: Cosmological Constraints. Astrophys. J. 2022, _938_ , 110.
  * [70] Sahni, V.; Saini, T.D.; Starobinsky, A.A.; Alam U. Statefinder-A new geometrical diagnostic of dark energy. JETP Lett. 2003, _77_ , 201.
  * [71] Sahni, V.; Shafieloo, A.; Starobinsky, A.A. Two new diagnostics of dark energy. Phys. Rev. D 2008, _78_ , 103502.
  * [72] Raychaudhuri, A. Relativistic Cosmology I.Phys. Rev. 1955, _98_ , 1123.


# Learning Interpretable Heuristics for WalkSAT

Yannet Interian1 Sara Bernardini2 1University of San Francisco
2Royal Holloway University of London
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Local search algorithms are well-known methods for solving large, hard
instances of the satisfiability problem (SAT). The performance of these
algorithms crucially depends on heuristics for setting noise parameters and
scoring variables. The optimal setting for these heuristics varies for
different instance distributions. In this paper, we present an approach for
learning effective variable scoring functions and noise parameters by using
reinforcement learning. We consider satisfiability problems from different
instance distributions and learn specialized heuristics for each of them. Our
experimental results show improvements with respect to both a WalkSAT baseline
and another local search learned heuristic.

##  1 Introduction

The satisfiability problem (SAT), one of the most studied NP-complete problems
in computer science, consists in determining if there exists an assignment
that satisfies a given Boolean formula. SAT algorithms typically assume that
formulas are expressed in conjunctive normal form (CNF). A CNF formula is a
conjunction of clauses; a clause is a disjunction of literals; and a literal
is a variable or its negation. SAT has a wide range of practical applications,
including electronic design automation, planning, scheduling and hardware
verification.

Stochastic local search (SLS) algorithms are well-known methods for solving
hard, large SAT instances (?). They are incomplete solvers: they typically run
with a pre-set number of iterations, after which they produce a valid
assignment or return “unsolved.” Algorithm 1 shows the pseudo-code of a
generic SLS algorithm. Like most SLS solvers, it starts by generating a random
assignment. If the formula is satisfied by this assignment, a solution is
found. Otherwise, a variable is chosen by a variable selection heuristic
(_pickVar_ in Algorithm 1) and that variable is flipped. The loop is repeated
until a solution is found or the maximum number of iterations is reached.

Algorithm 1 _SLS_ algorithm

Input: A formula $F$ in CNF form
Parameter: max_flips, max_tries
Output: If found, a satisfying assignment

for $i=1$ to max_tries do

$flips=0$

$X$ be a random initial assignment

while $flips$ $<$ max_flips do

if $X$ satisfies $F$ then

return $X$

Pick a variable $x$ using $pickVar$

$X$ $\leftarrow$ $flipVar(x)$

$flips$ ++

return unsolved

WalkSAT (?; ?; ?) and other successful local search algorithms select the
variable to flip from an unsatisfied clause (see Algorithm 2). After picking a
random unsatisfied clause $c$, the choice of which variable in $c$ to flip is
made in two possible ways: either a random variable is chosen, or a scoring
function is used to select the best variable to flip. The version of WalkSAT
in Algorithm 2 picks a variable with the smallest “break” value, where
$break(x)$ of a variable $x$ given an assignment $X$ is the number of clauses
that would become false by flipping $x$. Other algorithms and other versions
of WalkSAT use different heuristics (?; ?) for choosing $x$.

WalkSAT-type algorithms also use a noise parameter $p$ (see Algorithm 2) to
control the degree of greediness in the variable selection process. This
parameter has a crucial impact on the algorithms’ performance (?; ?; ?; ?).
Hoos et al. (?) propose a dynamic noise adaptation algorithm in which high
noise values are only used when the algorithms appear to not be making
progress.

Algorithm 2 _pickVar_ for WalkSAT

Pick a random unsatisfied clause $c$

if $rand()<p$ then

Pick a random variable $x$ from $c$

else

Pick a variable $x$ from $c$ with the smallest break value

return $x$

Designing SLS algorithms requires substantial problem-specific research and a
long trial-and-error process by the algorithm experts. Also, algorithms seldom
exploit the fact that real-world problems of the same type are solved again
and again on a regular basis, maintaining the same combinatorial structure,
but differing in the data. Problems of this type include, for example, SAT
encodings of AI Planning instances (?) and Bounded Model Checking instances
(?).

Recently, there has been increased interest in applying machine learning
techniques to design algorithms to tackle combinatorial optimization problems
(?; ?; ?; ?). In line with this work, our paper focuses on using machine
learning to design algorithms for SAT. More specifically, we investigate the
use of reinforcement learning to learn both adaptive noise strategies and
variable scoring functions for WalkSAT-type algorithms. We call the resulting
strategy _LearnWSAT_. The main contributions of this paper are as follows:

  * •

Our technique automatically learns a scoring function and an adaptive noise
strategy for WalkSAT-type algorithms.

  * •

Our scoring functions are simple and interpretable. When coded efficiently,
they would have a running time per iteration similar to WalkSAT.

  * •

Our approach outperforms both a WalkSAT baseline algorithm and a previously
published learned SLS-type algorithm (?).

  * •

Our technique uses a “warm-up” strategy designed to substantially decrease
training time.

  * •

Our algorithm, when trained on a specific distribution, generalizes well to
both unseen instances and larger instances of the same distribution.

We remark that our goal in this paper is to show how reinforcement learning
could be leveraged to make WalkSAT-type algorithms more efficient and their
design more practical; we do not aim to offer the fastest WalkSAT
implementation, which we leave as future work. 111The implementation can be
found here https://github.com/yanneta/learning_heuristics_sat

##  2 Related Work

The literature regarding SAT is vast. We focus here only on the following two
topics, which are the most pertinent to our contribution.

###  2.1 Machine Learning for SAT

Guo et al. (?) give an in-depth survey of machine learning for SAT. In their
classification, our work falls into the category described as “modifying local
search solvers with learning modules”. There are two other works (?; ?) that
fall into the same category.

Yolcu and Poczos (?) use reinforcement learning with graph neural networks to
learn an SLS algorithm. The graph neural network takes a factor graph
associated which the SAT formula and the current assignment to score each
variable. Scoring each variable at every iteration incurs a large overhead,
which leads the authors to run experiments only on small SAT instances. Our
work is similar to Yolcu and Poczos’s (?) in that we also use a model to score
variables. On the other hand, our approach differs from theirs in four ways.
Our scoring model is a linear function of a small set of features, which is
simple and interpretable. At every iteration, we only score variables from one
unsatisfied clause, which makes our model much more scalable and practical.
Our features are able to encode time dependencies (e.g. last time a variable
was flipped). We learn a separate noise strategy.

Zhang et al. (?) propose a system (NLocalSAT) for guiding the assignment
initialization of an SLS solver with a neural network. Their model feeds the
CNF formula into a Gated Graph neural network for feature extraction. The
neural network predicts an assignment for the SAT formula. The model is
trained to predict a satisfying assignment. The output of the neural network
is used to initialize SLS solvers. Whereas NLocalSAT modifies the
initialization of the SLS algorithm, our algorithm modifies its internal loop.
Those two improvements are potentially compatible.

Selsam et al. (2018) trained a message-passing neural network called NeuroSAT
to predict the satisfiability (SAT) or unsatisfiability (UNSAT) of problem
instances. The authors trained and evaluated NeuroSAT on random problem
instances that are similar to the ones used in our paper. NeuroSAT achieved an
accuracy of 85% and successfully solved 70% of the SAT problems. It is worth
noting that our approach focuses on predicting satisfiability and does not
directly address unsatisfiability. However, our approach demonstrates a
significantly higher accuracy on SAT instances.

###  2.2 Stochastic Local Search for SAT

Various strategies have been proposed for picking the variables to flip within
WalkSAT.

McAllester et al. (?) analyze six strategies. In all the strategies, a random
unsatisfied clause $c$ is selected, and the variable is chosen within $c$.
With probability $p$, a random variable is selected from $c$; otherwise, one
of the six following strategies is implemented. 1) Pick the variable that
minimizes the number of unsatisfiable clauses. 2) Pick the variable that
minimizes the break value (Algorithm 2). 3) Same as the previous strategy, but
never make a random move if one with break value 0 exits. 3) Pick the variable
that minimizes the number of unsatisfied clauses, but refuse to flip any
variable that has been flipped in the last $t$ steps. 5) Sort the variables by
the total number of unsatisfied clauses, then pick the one with the smallest
value. Break ties in favor of the least recently flipped variable. 6) Pick a
variable using a combination of least recently picked variable and number of
unsatisfied clauses.

ProbSAT (?) uses a scoring function based on the values $make(x)$ and
$break(x)$ and samples the variable to pick based on that scoring function.
Given a variable $x$ and an assignment $X$, $make(x)$ is the number of clauses
that would become true by flipping $x$. Note that $make(x)-break(x)$ is the
number of unsatisfiable clauses after flipping $x$. Balint and Schoning (?)
experiment with various types of scoring functions based on $make$ and $break$
and find that $make$ values can be ignored.

Hoos (?) proposes a dynamic noise strategy that uses higher values of noise
only when the algorithm is in an “stagnation” stage, which is when there is no
improvement in the objective function’s value over the last $\frac{m}{6}$
search steps, where $m$ is the number of clauses of the given problem
instance. Every incremental increase in the noise value is realized as
$p\leftarrow 0.8p+0.2$; the decrements are defined as $p\leftarrow 0.6p$ where
$p$ is the noise level.

The work by McAllester et al. (?) inspired our selection of features for the
variable ranking, and the paper by Balint and Schoning (?) led us to use
features based on $break(x)$ and ignore $make(x)$. Finally, the work in Hoos
(?) inspired us to learn an automated noise strategy.

##  3 Methodology

Algorithm 3 shows the pseudo-code for our $pickVar$ module. Our objective is
to learn the functions $p_{w}$ and $f_{\theta}$ in such a way that they
minimize the number of flips needed to solve a SAT problem. We now describe
these functions in detail.

Algorithm 3 _pickVar_ for LearnWSAT

Pick a random unsatisfied clause $c$

if $rand()<p_{w}$ then

Pick a random variable $x$ from $c$

else

Compute score $f_{\theta}(z)$ for each variable $z$ in $c$

$x\leftarrow$ sample $z$ with prob $\frac{e^{f_{\theta}(z)}}{\sum_{y\in
c}e^{f_{\theta}(y)}}$

return $x$

###  3.1 Variable Representation

To score each variable, we first compute some features that represent the
state of the variable at the current iteration $t$. From our discussion of
previous work in Section 2.2, we know that $break(x)$ is an important feature
in deciding the score of a variable. We also know, from previous work, that we
want to avoid flipping variables back and forth. We design features encoding
that information.

Let $age_{1}(x)$ be the last iteration in which $x$ was flipped and
$age_{2}(x)$ the last iteration in which $x$ was flipped and selected by the
algorithm using $f_{\theta}(x)$. Let $last_{K}(x)=1$ if $x$ was flipped in the
last $K$ iterations by $f_{\theta}(x)$. Let $\tilde{x}=min(x,10)$.

Based on this notation, we represent each variable via the following features:

  * •

$bk(x)=\log(1+break(\tilde{x}))$

  * •

$\Delta_{1}(x)=1-\frac{age_{1}(x)}{t}$

  * •

$\Delta_{2}(x)=1-\frac{age_{2}(x)}{t}$

  * •

$last_{5}(x)$

  * •

$last_{10}(x)$

We use $\tilde{x}$ and $\log$ in the feature $bk(x)$ to make the feature
independent of the size of the formulas. $bk(x)$ it is also normalized to be
between 0 and 1.

We have selected these features based on an extensive preliminary evaluation
performed on a variety of features and formulas. It would be easy to expand
our technique to include additional features whenever relevant.

Let
$\mathbf{f}(x)=(bk(x),\Delta_{1}(x),\Delta_{2}(x),last_{5}(x),last_{10}(x))$
be the vector representing the variable $x$ at iteration $t$, given a current
assignment $X$ for a formula $F$. Note that, to compute the vector, we keep
updating variables $age_{1},age_{2},last_{10}$, which is very cheap. Similar
to WalkSAT, $break(x)$ is only computed for variables on one clause at each
iteration.

###  3.2 Models for Scoring Variables and Controlling Noise

Our goal is to make our algorithm interpretable and fast, so we use a linear
model for scoring variables. Given a feature vector $\mathbf{f}=\mathbf{f}(x)$
for a variable $x$, $f_{\theta}(x)$ is a linear model on $\mathbf{f}$:

$f_{\theta}(x)=\theta_{0}+\sum_{i}\theta_{i}\cdot\mathbf{f}_{i}$

Inspired by the dynamic noise strategy discussed in Section 2.2, we define the
stagnation parameter $\delta$ as the number of iterations since the last
improvement in the number of satisfied clauses, divided by the number of
clauses. Instead of increasing or decreasing it at discrete intervals as in
Hoos (?), our noise is a continuous function of $\delta$, defined as

$p_{w}(\delta)=0.5\cdot Sigmoid(w_{0}+w_{1}\delta+w_{2}\delta^{2})$

We use the sigmoid function to ensure $p_{w}$ being between $0$ and $0.5$.
Those are commonly used values for noise. Parameters $w_{0},w_{1},w_{2}$ are
learned together with parameters $\\{\theta_{i}\\}_{i=0}^{5}$ by using
reinforcement learning.

After running our initial experiments, we noticed that the effect of the
stagnation parameter $\delta$ was almost negligible. Therefore, in most of our
experiments, we use a noise parameter that is a constant learned for each
instance distribution, that is, $p_{w}=0.5\cdot Sigmoid(w_{0})$.

###  3.3 Simplicity and Interpretability of Models

Domingos (?) states that one interpretation of Occam’s razor in machine
learning is the following: “ _Given two models with the same generalization
error, the simpler one should be preferred because simplicity is desirable in
itself._ ”

Following this basic principle, in our technique, we use simple functions
(linear and sigmoid functions) involving a small set of input variables and
show that we get better results than related algorithms that use much more
complex models, e.g. Yolcu and Poczos’s one (?). Simplicity is also valuable
because simple linear models are very fast to evaluate, which is crucial to
practical SAT solvers.

Interpretability refers to a model’s capacity to be “ _explained or presented
in understandable terms to a human_ ” (?). Linear models that use only a few
simple variables are typically considered highly interpretable. Our variable-
scoring model, which has just six coefficients, is therefore highly
interpretable. The interpretability of a model is useful because it allows us
to identify which features are significant and important and thus make
decisions about adding or subtracting features. If a feature has a coefficient
close to 0, we can infer that the feature lacks statistical significance and
should be eliminated. By providing insight into the impact of each model
feature, interpretability can help algorithm designers simplify the process of
adding, removing, and designing features.

Table 1 provides an example of the scoring parameters associated with random
3-SAT formulas of various sizes. The absolute value of each coefficient in the
table allows us to gauge the contribution of each variable to the model. As
demonstrated by the coefficients in Table 1, the $bk(x)$ feature has a notably
negative impact on the variable score, indicating its strong influence
compared to other features. Conversely, the coefficients associated with the
noise function $p_{w}(\delta)$ showed that $\delta$ was not a crucial feature,
allowing us to simplify our assumptions regarding the noise parameter. This
kind of insight can be extremely valuable.

size | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | $\theta_{4}$ | $\theta_{5}$ | $\theta_{0}$
---|---|---|---|---|---|---
$rand_{3}$ |  |  |  |  |  |
$(50,213)$ | -21.1 | -1.8 | -2.9 | -0.9 | -1.3 | 0.1
$(75,320)$ | -19.0 | -1.8 | -2.3 | -0.8 | -1.1 | 0.5
$(100,426)$ | -18.1 | -1.7 | -2.0 | -1.2 | -1.4 | 0.6
$(200,852)$ | -19.4 | -2.4 | -2.6 | -1.0 | -1.5 | -0.2
$rand_{4}$ |  |  |  |  |  |
$(30,292)$ | -20.2 | -1.2 | -3.2 | 0.9 | -2.5 | 0.28
$(50,487)$ | -14.3 | -1.0 | -1.4 | 0.7 | -2.1 | -0.31
Table 1: Coefficients of the scoring variable model learned for
$rand_{k}(n,m)$. The first column specifies the size of the formulas $(n,m)$.
$\theta_{0}$ is the model’s bias.

###  3.4 Training with Reinforcement Learning

To learn heuristics by using reinforcement learning (?), we formalize local
search for SAT as a Markov Decision Process (MDP). For clarity, we describe
the MDP assuming that the noise parameter is 0, that is, the algorithm always
picks a variable $x$ from a random unsatisfied clause $c$ using features
$\mathbf{f}(x)$.

For each problem distribution $D$, we have an MDP represented as a tuple
$(\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\gamma)$ where:

  * •

$\mathcal{S}$ is the set of possible states. The state encodes the information
needed at iteration $t$ to pick a variable to flip. In our setting, a state is
a tuple $(X,c,\\{\mathbf{f}(x)\\}_{x\in c},t)$, where $X$ is our current
assignment, $c$ is a clause unsatisfied by $X$, and $\\{\mathbf{f}(x)\\}_{x\in
c})$ are the set of features for all variables in $c$ and $t$ is the current
step. The formula $F$ uniformly sampled from $D$ is also part of the state,
but it is fixed through the episode. There are also two end states:
$end_{sat}$ and $end_{unsolved}$.

  * •

$\mathcal{A}$ is the set of actions. Given a state
$s=(X,c,\\{\mathbf{f}(x)\\}_{x\in c},t)$, the set of actions corresponds to
picking a variable to flip from the state’s clause $c$.

  * •

$\mathcal{P}$ is the transition probability function, defining the probability
of going from a state-action pair $(s,a)$ to the next state $s^{\prime}$. Let
$s=(X,c,\\{\mathbf{f}(x)\\}_{x\in c},t)$ be our current state, we pick a
variable $x$ in $c$ with probability $\frac{e^{f_{\theta}(x)}}{\sum_{y\in
c}e^{f_{\theta}(y)}}$, which gets us $X^{\prime}$, the assignment obtained
from $X$ by flipping variable $x$. If $X^{\prime}$ satisfies the formula $F$,
we move to the $end_{sat}$ state. If the max number of steps is reached and
$X^{\prime}$ does not satisfy $F$, we move to $end_{unsolved}$. Otherwise, we
move to $(X^{\prime},c^{\prime},\\{\mathbf{f}(x)\\}_{x\in c^{\prime}},t+1)$,
where $c^{\prime}$ is a random unsatisfied clause by the new assignment
$X^{\prime}$.

  * •

$\mathcal{R}(s)$ is the immediate reward after transitioning to state $s$.
$\mathcal{R}(end_{sat})=1$ and $0$ otherwise.

  * •

$\gamma\in(0,1)$ is the discount factor, which we set to less than 1 to
encourage finding solutions in fewer steps.

We reformulate the problem of learning informative heuristics for SAT into the
problem of finding an optimal policy $\pi$ for the MDP described above. We use
the well-known REINFORCE algorithm (?). Our policy $\pi(s)$ is determined by
the function $f_{\theta}(x)$ that we use to sample the variable to flip based
on the feature vector of each variable.

At each training iteration, we sample a batch of formulas from the
distribution $D$ and generate trajectories for each formula. We accumulate the
policy gradient estimates from all trajectories and perform a single update of
the parameters. Algorithm 4 shows the pseudo-code of the REINFORCE algorithm
for the case of constant noise and batch size of one.

Algorithm 4 _Reinforce_

Input: Training set (Train_ds) from a problem distribution $D$, policy
$\pi_{\theta}$, discount rate $\gamma$, learning rate $\alpha$

Initialize parameters of the policy $\pi_{\theta}$

Wam-up $\pi_{\theta}$ by fitting it to WalkSAT scoring function

for $i=1$ to Epocs do

Initialize $g\leftarrow 0$;

for $j=1$ to len(Train_ds) do

$F=$ Train_ds[j]

Init random assignment $X$; state $S_{0}$; $history=()$

for $t=0$ to max_flips do

if $X$ satisfies $F$ then

break

Sample action $a\sim\pi_{\theta}(S_{t})$

Append $(S_{t},a)$ to $history$

$X\leftarrow flip(X,a)$

Update state $S_{t+1}$

Set reward $r=1$ if $X$ satisfies $F$ and 0 otherwise.

for $t=0$ to $T=len(history)$ do

$(S_{t},a)\leftarrow history(t)$

$\hat{p}\leftarrow\pi_{\theta}(S_{t})$

$g\leftarrow g+\gamma^{T-t}r\nabla\log\hat{p}(a)$

$\theta\leftarrow\theta+\alpha g$

###  3.5 Training with a Warm-Up Strategy

By performing an extensive experimental evaluation, we found that the training
of our algorithm takes too long for formulas with over 50 variables when using
completely random heuristics and not initially finding a satisfying
assignment. Trials without satisfying assignments are not useful for training
since they have a reward of zero. To cope with this problem, we design a warm-
up strategy to speed up the training process. For a few epochs, we train the
function $f_{\theta}$ in such a way that the sampling mimics the pickVar
strategy from WalkSAT with probability $\frac{e^{f_{\theta}(z)}}{\sum_{y\in
c}e^{f_{\theta}(y)}}$. We cast this as a classification problem and use log-
loss and gradient descent to train $f_{\theta}$. Figure 1 displays the
training with and without warm-up for formulas in $rand_{3}(75,320)$, showing
the benefit of our approach.

Figure 1: Comparing median flips (log-scale) over epochs on training data for
LearnWSAT with and without a 5 epoch warm-up. Training with 1800 formulas of
$rand_{3}(75,320)$.

##  4 Experimental Setting

###  4.1 Data

We perform experiments using random formulas generated from the following
problems: random 3-SAT, random 4-SAT, clique detection, graph coloring and
dominating set. These distributions, except for random 4-SAT, are used in the
evaluation of GnnSLS by Yolcu and Poczos (?). To facilitate comparison, we use
the same problem distributions. They also used a vertex covering problem that
the CNFgen package (?) no longer supports, so we do not include this problem
in our experiments.

It has been observed empirically that random $K$-SAT problems are hard when
the problems are critically constrained, i.e. close to the SAT/UNSAT phase
boundary (?; ?). These problems are used as common benchmarks for SAT. The
threshold for $3$-SAT is when problems have roughly $4.26$ times as many
clauses as variables. To generate hard problems for random $4$-SAT, we set the
number of clauses to be $9.75$ times the number of variables (?). The other
three problems are NP-complete graph problems. For each of these problems, a
random Erdos–Rényi graph $G(N,p)$ is sampled. To sample from $G(N,p)$, a graph
with $N$ nodes is generated by sampling each edge with probability $p$.

For all these problem distributions, we generate random instances and keep
those that are satisfiable. We use the CNFgen package (?) to generate all
instances and Minisat (?) to filter out the unsatisfiable formulas.

###  4.2 Algorithms

For comparison, we use the SLS algorithm learned via reinforcement learning
developed by Yolcu and Poczos (?), which we call GnnSLS, and follow the same
experimental setup. We also consider one of the WalkSAT versions, as described
in Selman et al. (?). Again, we follow Yolcu and Poczos (?) in using this
particular WalkSAT version.

We wrote our algorithms in Python and PyTorch, which does not make them
competitive with state-of-the-art SAT solvers with respect to running time.
Indeed, our goal in this paper is to explore the power of reinforcement
learning for formulating effective SAT heuristics. To this aim, we offer a
prototype algorithm that proves the concept. Although we do not try here to
beat highly-optimized current SAT solvers, our results suggest that our
technique has the potential to compete with them if written efficiently.

For each problem distribution, we generate 2500 satisfiable formulas. From
these, 500 are used for testing, 1900 for training and 100 for validation.

As metrics, we use the median of the median number of flips, the average
number of flips and the percentage of instances solved.

###  4.3 Training with Reinforcement Learning

We train GnnSLS as described in Yolcu and Poczos (?)’s paper and use their
code from the related GitHub repository. The paper uses curriculum learning,
where training is performed on a sequence of problems of increasing
difficulty. For example, to train problems for $rand_{3}(50,213)$, the authors
start by first training on $rand_{3}(5,21)$, using the resulting model to
subsequently train on $rand_{3}(10,43)$, $rand_{3}(25,106)$ and
$rand_{3}(50,213)$.

As mentioned above, for experiments with random formulas, our models are
trained using 1900 instances. The 100 validation instances are used to select
the model with the best median number of steps. We train for 60 epochs using
one cycle training (?) and AdamW (?) as the optimizer (a link to our GitHub
repository will be provided in due course). Most of our experiments are run
with a discount factor of 0.5.

###  4.4 Evaluation

For evaluation, we use $max\\_tries=10$ and $max\\_flips=10000$ unless
otherwise specified. As said above, for randomly generated problems, we use
500 instances for testing. The noise probability for WalkSAT and GnnSLS is set
to $p=\frac{1}{2}$ as in the experiments by Yolcu and Poczos (?).

| LearnWSAT | GnnSLS | WalkSAT
---|---|---|---
$rand_{3}(50,213)$ |  |  |
m-flips | 119 | 352 | 356
a-flips | 384 | 985 | 744
solved | 100% | 99.6% | 100%
$color_{5}(20,0.5)$ |  |  |
m-flips | 103 | 137 | 442
a-flips | 225 | 497 | 787
solved | 100% | 100% | 100%
$clique_{3}(20,0.05)$ |  |  |
m-flips | 68 | 200 | 176
a-flips | 91 | 345 | 238
solved | 100% | 100% | 100%
$domeset_{4}(12,0.2)$ |  |  |
m-flips | 65 | 72 | 171
a-flips | 97 | 242 | 288
solved | 100% | 100% | 100%
$rand_{4}(50,487)$ |  |  |
m-flips | 685 | - | 2044
a-flips | 1484 | - | 3302
solved | 100% | - | 96%
Table 2: Performance of LearnWSAT compared to GnnSLS and WalkSAT. Three
metrics are presented: median (m-flips) and average number of flips (a-flips),
and percentage solved (solved).

##  5 Experimental Results

#### Comparison to GnnSLS and WalkSAT.

Table 2 summarizes the performance of LearnWSAT compared to GnnSLS and
WalkSAT. We present results for five classes of problems, $rand_{3}(50,213)$,
$rand_{4}(30,292)$, $color_{5}(20,0.5)$, $clique_{3}(20,0.05)$ and
$domeset_{4}(12,0.2)$ and three metrics, median number of flips (m-flips),
average number of flips (a-flips), and percentage solved (solved). Table 3
indicates the number of variables and clauses in the sampled formulas and
gives a sense of the size of the SAT problems we tackle. Table 2 shows that,
after training, LearnWSAT requires substantially fewer steps than GnnSLS and
WalkSAT to solve the respective problems.

Distribution | $n$ | $m$
---|---|---
$rand_{k}(n,m)$ | $n$ | $m$
$color_{5}(20,0.5)$ | 100 | 770
$clique_{3}(20,0.05)$ | 60 | 1758
$domeset_{4}(12,0.2)$ | 60 | 996
Table 3: Size of the formula used in our evaluation. The distribution
$rand_{k}(n,m)$ has exactly $n$ variables and $m$ clauses. For all the other
distributions, we show the maximum number of variables $n$ and clauses $m$ in
the sampled formulas.

Our algorithm performs better than WalkSAT because it optimizes the variable
scoring and the noise parameter to the particular distribution of SAT
problems. Our technique is also better than GnnSLS because of the following
two reasons. First, we speculate that GnnSLS underfits the problem. The SAT
encoding and the model used by GnnSLS are more sophisticated but also much
more complex than our approach. It is not possible to directly train the
GnnSLS algorithm with problems that have a few variables (e.g. 50 variables).
To get the GnnSLS encoding to work well, smarter training and more data are
needed. Second, our approach uses time-dependent variables (the last time a
variable has been flipped), which GnnSLS is unable to encode.

#### Generalization to larger instances.

In Table 4, we compare the performance of LearnWSAT trained on data sets of
different sizes to assess how well the algorithm generalizes to larger
instances after having been trained on smaller ones. We consider random 3-SAT
instances of different sizes, $rand_{3}(n,m)$. As in Table 2, we consider
three metrics: median number of flips (m-flips), average number of flips
(a-flips), and percentage solved (solved). The second column reports the
performance of LearnWSAT (indicated LWSAT for brevity) on instances of
different sizes when the algorithm is trained on $rand_{3}(50,213)$ only. In
the third column, for comparison, we report the performance of LearnWSAT when
it is trained and evaluated on instances of the same size. The fourth column
reports the performance of GnnSLS when the algorithm is trained on
$rand_{3}(50,213)$ only. Finally, the last column reports the WalkSAT
(indicated WSAT) baseline.

The table shows that our model evaluated on $rand_{3}(50,213)$ performs
similarly or better than models trained on larger instances. Training becomes
much more expensive as a function of the size of the formula, but this result
suggests that we can train on smaller formulas of the same distribution.
GnnSLS trained on smaller instances can also be evaluated on larger problems
of the same distribution, but the results seem to degrade as the formulas get
larger.

| LWSAT | LWSAT | GnnSLS | WSAT
---|---|---|---|---
| (50, 213) |  | (50, 213) |
$(50,213)$ |  |  |  |
m-flips | 119 | 119 | 352 | 356
a-flips | 384 | 384 | 985 | 744
solved | 100% | 100% | 99.6% | 100%
$(75,320)$ |  |  |  |
m-flips | 260 | 286 | 969 | 880
a-flips | 904 | 948 | 2253 | 1772
solved | 100% | 100% | 96.6% | 98%
$(100,426)$ |  |  |  |
m-flips | 503 | 575 | 2264 | 1814
a-flips | 1650 | 1682 | 3816 | 3132
solved | 100% | 100% | 85.6% | 93%
$(200,852)$ |  |  |  |
m-flips | 4272 | 4005 | 10000 | 10000
a-flips | 5329 | 5085 | 8359 | 7497
solved | 96.2% | 95.6% | 26% | 46.2%
Table 4: Performance of our algorithm evaluated on different instances of the
same distribution. We consider $rand_{3}(n,m)$ formulas of different sizes
($n$ and $m$ refer to the number of variables and clauses in the sampled
formulas). The second column corresponds to evaluating our algorithm
(indicated as LWSAT) trained on formulas from $rand_{3}(50,213)$ only. The
third column corresponds to evaluating the algorithm on instances of the same
size used for training. The fourth corresponds to evaluating GnnSLS on the
algorithm trained on $rand_{3}(50,213)$. The last column reports the WalkSAT
(indicated WSAT) baseline. We consider three metrics: median (m-flips) and
average number of flips (a-flips), and percentage solved (solved).

Table 5 shows results on instances that are harder than the ones shown before.
In particular, Minsat is not able to solve some of the instances of
$rand_{3}(500,2130)$ and $rand_{4}(200,1950)$ in less than ten hours. We
generated 100 problems from $rand_{3}(300,1278)$, $rand_{3}(500,2130)$ and
$rand_{4}(200,1950)$, respectively.

These instances are generated at the SAT/UNSAT threshold, therefore around 50%
of them are supposed to be satisfiable. In the case of $rand_{4}(200,1950)$,
it seems that a few more are satisfiable since LearnWSAT is able to solve 68%
of them.

| LearnWSAT | WalkSAT
---|---|---
$rand_{3}(300,1278)$ | 48% | 26%
$rand_{3}(500,2130)$ | 36% | 9%
$rand_{4}(200,1950)$ | 68% | 0%
Table 5: Performance of our algorithm evaluated on larger instances of
$rand_{k}(n,m)$. These instances have not been checked for satisfiability.
Around 50% of them are expected to be satisfiable. The metric shown is
percentage solved. The max_flip parameter is set to 50000 for both solvers.
LearnWSAT was trained on $rand_{3}(50,213)$ for the first two rows and on
$rand_{4}(50,487)$ for the last row.

#### Noise parameter.

In our initial experiments, we learned a noise function that depended on the
stagnation parameter $\delta$. After inspecting the function, we noticed that
the effect of $\delta$ is negligible. In Figure 2, we show the learned noise
function as used by the algorithm at evaluation time.

Figure 2: The data corresponding to each line comes from running an evaluation
on a single formula for each of the five problem distributions
($color_{5}(20,0.5)$, $clique_{3}(20,0.05)$, $domeset_{4}(12,0.2)$,
$rand_{3}(50,213)$, $rand_{4}(30,292)$) until the SAT assignment is found. The
plot shows the noise parameter as a function of the iteration.

We plot the noise function against the iteration until the formula is solved.
The stagnation parameter varies per iteration, but these curves show very
little variation. We ran experiments in which we fixed $p_{w}$ to be a
constant dependent on the distribution and found that the results are similar
to when the noise function depends on $\delta$. In particular, we optimize
$p_{w}=0.5\cdot Sigmoid(w)$ by finding a single parameter $w$ per
distribution. After these initial experiments, we ran all the others (as they
are reported here) with fixed constants. Note that these constants are small
compared to typical values used for WalkSAT ($p=1/2$). This is because our
$PickVar$ algorithm (shown in Algorithm 3) injects noise by sampling instead
of deterministically picking variables as in the original $PickVar$ algorithm
of WalkSAT (Algorithm 2).

#### Impact of the discount factor.

We ran experiments to understand the dependencies of our results on the value
of the discount factor for reinforcement learning. Figure 3 shows the median
flips as a function of the discount factor. The gray area shows the confidence
intervals for each curve. We find that various discount factors give similar
results.

Figure 3: Comparing median flips as a function of the discount factor for
various datasets. The lines correspond to training and evaluation on instances
of the following distributions: $rand_{3}(50,213)$, $rand_{4}(30,292)$,
$color_{5}(20,0.5)$, $clique_{3}(20,0.05)$ and $domeset_{4}(12,0.2)$

#### Impact of the size of training data.

Figure 4 shows the median flips as a function of the size of the training
data. The experiment uses formulas from $rand_{3}(50,213)$. The plot shows
that we need a training size of at least 40 to learn an algorithm that is
better than WalkSAT. For optimal results, we need at least 160 formulas. To
run the experiments with smaller datasets, we increased the number of warm-up
steps from 5 to 50 and the amount of epochs from 60 to 200.

Figure 4: Comparing median flips as a function of the size of the training
data. Training with $rand_{3}(50,213)$ formulas. Gray region corresponds to
confidence intervals. Gray line corresponds to median flips for the WalkSAT
baseline.

##  6 Conclusions and Future Work

In this paper, we present LearnWSAT, a technique that discovers effective
noise parameters and scoring variable functions for WalkSAT-type algorithms.
Thanks to them, LearnWSAT uses substantially fewer flips than a WalkSAT
baseline, as well as an existing learned SLS-type algorithm, to solve the
satisfiability problem. Although we do not focus on optimizing the
implementation of LearnWSAT in this paper, our experiments suggest that, when
coded efficiently, our technique could compete with state-of-the-art solvers.

Despite improving over algorithms in the literature, we note that a limitation
of LearnWSAT is the need to pre-define a set of features. In addition,
training is slow for formulas with 150 variables or more. The last limitation
is mitigated by the fact that, as we have shown in the experiments, models
trained on smaller formulas generalize well to larger ones. Overcoming these
limitations is part of our future work.

Finally, we remark that the ideas presented in this work are general and could
be adapted to solve other hard combinatorial problems.

## References

  * 2012 Balint, A., and Schoning, U.  2012\.  Choosing probability distributions for stochastic local search and the role of make versus break.  In Cimatti, A., and Sebastiani, R., eds., Theory and Applications of Satisfiability Testing – SAT 2012.  Berlin, Heidelberg: Springer Berlin Heidelberg.
  * 2016 Bello, I.; Pham, H.; Le, Q. V.; Norouzi, M.; and Bengio, S.  2016\.  Neural combinatorial optimization with reinforcement learning.  arXiv preprint arXiv:1611.09940.
  * 2005 Benedetti, M., and Bernardini, S.  2005\.  Incremental compilation-to-sat procedures.  In Hoos, H. H., and Mitchell, D. G., eds., Theory and Applications of Satisfiability Testing, 46–58.  Berlin, Heidelberg: Springer Berlin Heidelberg.
  * 2021 Bengio, Y.; Lodi, A.; and Prouvost, A.  2021\.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  European Journal of Operational Research 290(2):405–421.
  * 1999 Domingos, P.  1999\.  The role of occam’s razor in knowledge discovery.  Data mining and knowledge discovery 3(4):409–425.
  * 2017 Doshi-Velez, F., and Kim, B.  2017\.  Towards a rigorous science of interpretable machine learning.  arXiv preprint arXiv:1702.08608.
  * 2003 Eén, N., and Sörensson, N.  2003\.  An extensible sat-solver.  In Giunchiglia, E., and Tacchella, A., eds., SAT, volume 2919 of Lecture Notes in Computer Science, 502–518.  Springer.
  * 1994 Gent, I. P., and Walsh, T.  1994\.  The sat phase transition.  In ECAI, volume 94, 105–109.  PITMAN.
  * 2022 Guo, W.; Yan, J.; Zhen, H.-L.; Li, X.; jie Yuan, M.; and Jin, Y.  2022\.  Machine learning methods in solving the boolean satisfiability problem.  ArXiv abs/2203.04755.
  * 1999 Hoos, H. H.  1999\.  On the run-time behaviour of stochastic local search algorithms for SAT.  In Hendler, J., and Subramanian, D., eds., Proceedings of the Sixteenth National Conference on Artificial Intelligence and Eleventh Conference on Innovative Applications of Artificial Intelligence, July 18-22, 1999, Orlando, Florida, USA, 661–666.  AAAI Press / The MIT Press.
  * 2002 Hoos, H. H.  2002\.  An adaptive noise mechanism for walksat.  In AAAI/IAAI.
  * 2009 Kautz, H. A.; Sabharwal, A.; and Selman, B.  2009\.  Incomplete algorithms.  In Biere, A.; Heule, M.; van Maaren, H.; and Walsh, T., eds., Handbook of Satisfiability, volume 185 of Frontiers in Artificial Intelligence and Applications. IOS Press.  185–203.
  * 2017 Khalil, E.; Dai, H.; Zhang, Y.; Dilkina, B.; and Song, L.  2017\.  Learning combinatorial optimization algorithms over graphs.  Advances in neural information processing systems 30\.
  * 2017 Lauria, M.; Elffers, J.; Nordström, J.; and Vinyals, M.  2017\.  Cnfgen: A generator of crafted benchmarks.  In Theory and Applications of Satisfiability Testing – SAT 2017, Lecture Notes in Computer Science, 464–473.  Germany: Springer.  20th International Conference on Theory and Applications of Satisfiability Testing, SAT 2017.
  * 2017 Loshchilov, I., and Hutter, F.  2017\.  Fixing weight decay regularization in adam.  CoRR abs/1711.05101.
  * 1997 McAllester, D. A.; Selman, B.; and Kautz, H. A.  1997\.  Evidence for invariants in local search.  In Kuipers, B., and Webber, B. L., eds., Proceedings of the Fourteenth National Conference on Artificial Intelligence and Ninth Innovative Applications of Artificial Intelligence Conference, AAAI 97, IAAI 97, July 27-31, 1997, Providence, Rhode Island, USA, 321–326.  AAAI Press / The MIT Press.
  * 1992 Mitchell, D.; Selman, B.; Levesque, H.; et al.  1992\.  Hard and easy distributions of sat problems.  In Aaai, volume 92, 459–465.
  * 2008 Robinson, N.; Gretton, C.; Pham, D. N.; and Sattar, A.  2008\.  A compact and efficient sat encoding for planning.  In ICAPS, 296–303.
  * 1993 Selman, B.; Kautz, H. A.; and Cohen, B.  1993\.  Local search strategies for satisfiability testing.  In Johnson, D. S., and Trick, M. A., eds., Cliques, Coloring, and Satisfiability, Proceedings of a DIMACS Workshop, New Brunswick, New Jersey, USA, October 11-13, 1993, volume 26 of DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 521–531.  DIMACS/AMS.
  * 1994 Selman, B.; Kautz, H. A.; and Cohen, B.  1994\.  Noise strategies for improving local search.  In Hayes-Roth, B., and Korf, R. E., eds., Proceedings of the 12th National Conference on Artificial Intelligence, Seattle, WA, USA, July 31 - August 4, 1994, Volume 1, 337–343.  AAAI Press / The MIT Press.
  * 1996 Selman, B.; Mitchell, D. G.; and Levesque, H. J.  1996\.  Generating hard satisfiability problems.  Artificial intelligence 81(1-2):17–29.
  * 2017 Smith, L. N., and Topin, N.  2017\.  Super-convergence: Very fast training of residual networks using large learning rates.  CoRR abs/1708.07120.
  * 2018 Sutton, R. S., and Barto, A. G.  2018\.  Reinforcement Learning: An Introduction.  The MIT Press, second edition.
  * 1992 Williams, R. J.  1992\.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  Machine learning 8(3):229–256.
  * 2019 Yolcu, E., and Póczos, B.  2019\.  Learning local search heuristics for boolean satisfiability.  Advances in Neural Information Processing Systems 32\.
  * 2020 Zhang, W.; Sun, Z.; Zhu, Q.; Li, G.; Cai, S.; Xiong, Y.; and Zhang, L.  2020\.  Nlocalsat: Boosting local search with solution prediction.  CoRR abs/2001.09398.