Accelerated Fourier SAT (AFSAT): Fully Realising a GPU-based Symmetric Pseudo-Boolean SAT Solver

# Accelerated Fourier SAT: Fully Realising a GPU-based Symmetric Pseudo-Boolean SAT Solver
Source: [https://arxiv.org/html/2606.06641](https://arxiv.org/html/2606.06641)
\\hideLIPIcs

School of Computing, Australian National University, Canberra, Australia\.cody\.christopher@anu\.edu\.auhttps://orcid\.org/0000\-0001\-8444\-2292 School of Computing, Australian National University, Canberra, Australia\.charles\.gretton@anu\.edu\.auhttps://orcid\.org/0000\-0001\-9803\-0168\\CopyrightCody J Christopher and Charles Gretton\\ccsdesc\[500\]Mathematics of computing Combinatorial optimization\\ccsdesc\[500\]Mathematics of computing Solvers\\ccsdesc\[500\]Mathematics of computing Nonconvex optimization\\supplement\\afsimplementation withApache\-2\.0/GPL\-2\.0\-or\-laterlicenses:\\supplementdetails\[linktext=, cite=, subcategory=Source Code, swhid= \]Softwarehttps://github\.com/cjchristopher/accelerated\-fourier\-sat\\EventEditors\\EventNoEds0\\EventLongTitle\\EventShortTitle\\EventAcronym\\EventYear\\EventDate\\EventLocation\\EventLogo\\SeriesVolume\\ArticleNo

###### Abstract

We present Accelerated Fourier SAT \(\\afs\), a GPU\-accelerated solver for pseudo\-Boolean satisfiability based on continuous local search \(CLS\)\.\\afsrealises the proof\-of\-concept approach,FastFourierSAT, into a fully\-engineered solver supporting any heterogeneous mixture of symmetric constraint types and lengths within a single problem instance\. Using theJAXcompiler,\\afsleverages pure function composition, automatic vectorisation, automatic differentiation, and just\-in\-time \(JIT\) compilation to perform massively parallel CLS across batches of candidate assignments\. We demonstrate substantially improved numerical stability, runtime performance, and memory efficiency over the proof\-of\-concept\. We achieve this by way of identifying and addressing various limitations that arise from memory latency and floating\-point representation, as well as leveraging automatic parallelisation and compact representations\. The inherent representational and stability limitations of floating point are partially addressed by a tailored discrete Fourier transform implementation\. We achieve near\-linear throughput when scaling to multiple accelerators viaJAXarray sharding\.

###### keywords:

Satisfiability, pseudo\-Boolean, SAT Solver, continuous local search, combinatorial optimization, hardware acceleration

## 1Introduction

Continuous local search \(CLS\) offers a compelling search paradigm for solving satisfiability \(SAT\) problems that are expressed using symmetric pseudo\-Boolean \(PB\) constraints\. The approach relaxes Boolean problem variables to real\-valued variables by way of the Walsh\-Fourier transform from Boolean function analysis\[odonnell14\]\. Revisiting this approach, SAT is reformulated as a bounded continuous optimisation problem amenable to gradient\-based search methods\. This approach was formalised for theFourierSATproof\-of\-concept\[kyrillidis2021solving\]and subsequently extended for parallelised GPU computation inFastFourierSAT\[cen2025massively\], which demonstrated that the Walsh\-Fourier expansion can be evaluated efficiently using a vectorised Discrete Fourier Transform \(DFT\) inO∗​\(log⁡k\)O^\{\*\}\(\\log k\)time wherekkis the number of variables \(typically also literals\) in a constraint \(clause\), andO∗​\(⋅\)=d​e​fOp→∞​\(Tp​\(⋅\)\)O^\{\*\}\(\\cdot\)\\stackrel\{\{\\scriptstyle def\}\}\{\{=\}\}O\_\{p\\rightarrow\\infty\}\(T\_\{p\}\(\\cdot\)\)is idealised parallel execution time with infinite resources\.

We present Accelerated Fourier SAT \(\\afs\), a ground\-up re\-engineered and extended implementation of this CLS approach\.\\afsprovides the following improvements and novel contributions as a tool:

- •Heterogeneous constraint support\.\\afsis the first CLS implementation to support problems with an arbitrary mix of common PB constraint types of varying lengths within a single problem instance\.
- •Improved performance and efficiency\.We demonstrate better execution times \(up to parity in the worst case\) with less overhead, substantially reduced GPU memory consumption, and higher peak throughput \(evaluations/searches per unit time\) compared toFastFourierSAT\.
- •Multi\-GPU scaling\.We demonstrate near\-linear scaling across multiple GPUs using distributed array sharding inJAX\[deepmind2020jax\], in a single\-program\-multiple\-data \(SPMD\), or compute\-follows\-data paradigm\.
- •Numerical stability and precision improvements\.We identify and address sources of floating\-point errors and representational deviations through a tailored DFT matrix construction with deferred division\. This establishes a practical maximum constraint length of approximately 50 variables\.
- •Partial assignment integration\.\\afsaccepts partial variable assignments as input, opening up avenues to use it as a sub\-solver within decomposition\-based architectures such asDagster\[10\.1007/978\-3\-031\-20862\-1\_6\], or other portfolio approaches\.

## 2Background

### 2\.1Continuous Relaxation via Walsh\-Fourier Expansion

CLS operates on a continuous relaxation of a SAT problem\. For a Boolean formulaϕ\\phiovernnvariables, we seek a relaxation as a polynomial on the Boolean hypercube𝒬n\\mathcal\{Q\}^\{n\}that preserves satisfiability \(i\.e\. a solution to the relaxed formulation provides a solution for the discrete problem\)\. The \(Walsh\-\)Fourier expansion \(FE\)\[odonnell14\]of Boolean functions provides this relaxation, mapping\{True,False\}\\left\\\{\\texttt\{True\},\\texttt\{False\}\\right\\\}to\{−1,1\}\\left\\\{\-1,1\\right\\\}:

FEϕ​\(𝐗\)=def∑S∈2𝐗f^​\(S\)​∏xi∈Sxi\\texttt\{FE\}\_\{\\phi\}\(\\mathbf\{X\}\)\\stackrel\{\{\\scriptstyle\\text\{def\}\}\}\{\{=\}\}\\sum\_\{S\\in 2^\{\\mathbf\{X\}\}\}\\hat\{f\}\(S\)\\prod\_\{x\_\{i\}\\in S\}x\_\{i\}\(1\)wheref^​\(S\)\\hat\{f\}\(S\)are called the Fourier coefficients\. For*symmetric*constraints—those whose truth value depends only on the*simple count*of true literals—closed\-form solutions for the coefficients exist that are computable in polynomial time at worst\[kyrillidis2021solving\]\. Symmetric constraints necessarily have idempotent variable weights\.

#### 2\.1\.1Satisfaction as Optimisation

Given a formulaϕ\\phiinnnvariables, decomposed intommsymmetric constraintsC1,…,CmC\_\{1\},\\ldots,C\_\{m\}, the satisfiability problem is expressed as the bounded minimisation:

min𝐗​∑k=1mFECk​\(𝐗\)​subject to​𝐗∈𝒬n\\min\_\{\\mathbf\{X\}\}\\sum\_\{k=1\}^\{m\}\\texttt\{FE\}\_\{C\_\{k\}\}\(\\mathbf\{X\}\)\\text\{ subject to \}\\mathbf\{X\}\\in\\mathcal\{Q\}^\{n\}\(2\)where an assignment𝐗∈\[−1,1\]n\\mathbf\{X\}\\in\[\-1,1\]^\{n\}*satisfies*ϕ\\phiif∑kFECk​\(𝐗\)=−m\\sum\_\{k\}\\texttt\{FE\}\_\{C\_\{k\}\}\(\\mathbf\{X\}\)=\-m\. This formulation is differentiable, non\-convex, bounded, and saddle\-dense\. This motivates the use of projected and/or bounded gradient search methods as the first choice of search algorithms\. Observing thatFEis differentiable and fast parallelised automatic differentiation is available, we would also consider second\- or higher\-order methods, should they exist\.

#### 2\.1\.2DFT\-Based Vectorised Evaluation

A key insight from Cen et al\.\[cen2025massively\]is that the evaluation of the Fourier expansion can be performed using a DFT\. We observe that the terms inFEϕ\\texttt\{FE\}\_\{\\phi\}sharing a coefficientf^​\(S\)\\hat\{f\}\(S\)are all the size\|S\|\\left\\lvert S\\right\\rvertcombinations of the variables of the clause\. These particular sums are well\-studied polynomials known as the*elementary symmetric polynomials*\(ESPs\)\. Since the sequence of all lengthkkESPs up tonnvariables,𝒆kn\\bm\{e\}^\{n\}\_\{k\}, can be computed with \(0\-padded\) linear convolutions, we can compute the evaluation of a Fourier expanded formula with vectorised operations in the frequency domain by making use of the linear convolution theorem:

f^ϕ\\displaystyle\\hat\{f\}\_\{\\phi\}≡\[f^ϕ​\(∅\),f^ϕ​\(1\),…,f^ϕ​\(n\)\]\\displaystyle\\equiv\[\\hat\{f\}\_\{\\phi\}\(\\emptyset\),\\hat\{f\}\_\{\\phi\}\(1\),\\ldots,\\hat\{f\}\_\{\\phi\}\(n\)\]𝒆n\\displaystyle\\bm\{e\}^\{n\}≡\[e0n,e1n,…,enn\]=\(\[x1,1,0,…\]∗\[x2,1,0,…\]∗…∗\[xn,1,0,…\]\)\\displaystyle\\equiv\\left\[e\_\{0\}^\{n\},e\_\{1\}^\{n\},\\ldots,e\_\{n\}^\{n\}\\right\]=\\left\(\[x\_\{1\},1,0,\\ldots\]\*\[x\_\{2\},1,0,\\ldots\]\*\\ldots\*\[x\_\{n\},1,0,\\ldots\]\\right\)=W∗​W​\(\[x1,1,0,…\]∗\[x2,1,0,…\]∗…∗\[xn,1,0,…\]\)\\displaystyle=W^\{\*\}W\\left\(\[x\_\{1\},1,0,\\ldots\]\*\[x\_\{2\},1,0,\\ldots\]\*\\ldots\*\[x\_\{n\},1,0,\\ldots\]\\right\)=W∗​\(W​\(\[x1,1,0,…\]\)⋅W​\(\[x2,1,0,…\]\)⋅…⋅W​\(\[xn,1,0,…\]\)\)\\displaystyle=W^\{\*\}\\left\(W\\left\(\[x\_\{1\},1,0,\\ldots\]\\right\)\\cdot W\\left\(\[x\_\{2\},1,0,\\ldots\]\\right\)\\cdot\\ldots\\cdot W\\left\(\[x\_\{n\},1,0,\\ldots\]\\right\)\\right\)=W∗​\(\[1\+x1,ω\+x1,…,ωn\+x1\]⋅…⋅\[1\+xn,…,ωn\+xn\]\)\\displaystyle=W^\{\*\}\\left\(\[1\+x\_\{1\},\\omega\+x\_\{1\},\\ldots,\\omega^\{n\}\+x\_\{1\}\]\\cdot\\ldots\\cdot\[1\+x\_\{n\},\\ldots,\\omega^\{n\}\+x\_\{n\}\]\\right\)=W∗​\(\[1,ω,ω2,…,ωn\]T\+\[x1,x2,…,xn\]\)\\displaystyle=W^\{\*\}\\left\(\\left\[1,\\omega,\\omega^\{2\},\\ldots,\\omega^\{n\}\\right\]^\{T\}\+\\;\\left\[x\_\{1\},x\_\{2\},\\ldots,x\_\{n\}\\right\]\\right\)⊳\(Outer addition\)\\displaystyle\\mathllap\{\\rhd\\text\{\\small\(Outer addition\)\}\}=W∗​\(\[∏i=1n\(1\+xi\),∏i=1n\(ω\+xi\),…,∏i=1n\(ωn\+xi\)\]\)\\displaystyle=W^\{\*\}\\left\(\\left\[\\prod\_\{i=1\}^\{n\}\{\(1\+x\_\{i\}\)\},\\prod\_\{i=1\}^\{n\}\{\(\\omega\+x\_\{i\}\)\},\\ldots,\\prod\_\{i=1\}^\{n\}\{\(\\omega^\{n\}\+x\_\{i\}\)\}\\right\]\\right\)=W∗​\(\[∏i=1n\(ωj\+xi\)\]j=0n\)\\displaystyle=W^\{\*\}\\left\(\\left\[\\prod\_\{i=1\}^\{n\}\{\(\\omega^\{j\}\+x\_\{i\}\)\}\\right\]\_\{j=0\}^\{n\}\\right\)FEϕ\\displaystyle\\texttt\{FE\}\_\{\\phi\}≡∑\(f^ϕ⋅𝒆n\)\\displaystyle\\equiv\\sum\{\\left\(\\hat\{f\}\_\{\\phi\}\\cdot\\bm\{e\}^\{n\}\\right\)\}=∑\(f^ϕ​W∗⋅\(\[∏i=1n\(ωj\+xi\)\]j=0n\)\)\\displaystyle=\\sum\{\\left\(\\hat\{f\}\_\{\\phi\}W^\{\*\}\\cdot\\left\(\\left\[\\prod\_\{i=1\}^\{n\}\{\(\\omega^\{j\}\+x\_\{i\}\)\}\\right\]\_\{j=0\}^\{n\}\\right\)\\right\)\}Where∗\*is the convolution operator on sequences,⋅\\cdotis the Hadamard product \(element\-wise multiplication\),WWis the DFT matrix inn\+1n\+1dimensions \(W∗W^\{\*\}the conjugate transpose\), andω\\omegathe primitive\(n\+1\)th\(n\+1\)^\{\\text\{th\}\}root of unity\. By taking the DFT of the ESP sequence convolution, the computation can be expressed as element\-wise products and sums\.

## 3System Design and Implementation

### 3\.1Architecture

\\afs

is implemented in Python using theJAX\[deepmind2020jax\]and XLA compilers\. The architecture comprises:

1. \(a\)Problem ingestion: The specification and parsing of PB\-encoded problem instances specified in standard DIMACS or our hybrid variant \(Appendix[A](https://arxiv.org/html/2606.06641#A1)\)\.
2. \(b\)Fourier coefficient computation: Closed\-form computation of coefficients for all supported symmetric constraint types \(§[3\.2](https://arxiv.org/html/2606.06641#S3.SS2)\)\.
3. \(c\)DFT precision: Tailored calculation of roots of unity for conjugate symmetry \(§[3\.5\.1](https://arxiv.org/html/2606.06641#S3.SS5.SSS1)\)\.
4. \(d\)Compiled solver kernel: JIT\-compiled search algorithm with automatic differentiation, vectorised over batches of candidate assignments\.
5. \(e\)Multi\-GPU distribution: Sharding of batched assignments across available GPUs\.

The solver kernel is compiled once and executed repeatedly with different random initialisations\. Compilation produces an optimised XLA HLO program targeting the available accelerator architecture\.\\afsexposes many additional parameters and heuristic options which can be adjusted to target specific problem types\.

### 3\.2Supported Constraint Types

Unlike prior CLS implementations which are restricted to either a single constraint type, or a particular fixed set of constraint types per problem,\\afssupports heterogeneous problems containing any combination of the following symmetric pseudo\-Boolean constraint types:

Constraint TypePB FormCoefficient CostDisjunction \(OR\)∑nxn≥1\\sum\_\{n\}x\_\{n\}\\geq 1O​\(1\)O\(1\)At most one \(AMO\)∑nxn≤1\\sum\_\{n\}x\_\{n\}\\leq 1O​\(n\)O\(n\)Exactly one \(EO\)∑nxn=1\\sum\_\{n\}x\_\{n\}=1O​\(n\)O\(n\)Exactlykk\(EK\)∑nxn=k\\sum\_\{n\}x\_\{n\}=kO​\(n​log2⁡n\)O\(n\\log^\{2\}\{n\}\)Not all equal \(NAE\)⋀\{∑nxn<n∑nxn\>0\\bigwedge\\begin\{cases\}\\textstyle\\sum\_\{n\}x\_\{n\}<n\\\\ \\textstyle\\sum\_\{n\}x\_\{n\}\>0\\end\{cases\}O​\(1\)O\(1\)Exclusive Or \(XOR\)∑nxn≡1\(mod2\)\\sum\_\{n\}x\_\{n\}\\equiv 1\\pmod\{2\}O​\(n\)O\(n\)Cardinality\-kk\(CARD\)∑nxn≥k\\sum\_\{n\}x\_\{n\}\\geq kO​\(n​log2⁡n\)O\(n\\log^\{2\}\{n\}\)Table 1:Supported symmetric PB constraint types and asymptotic cost of computing their Walsh\-Fourier expansion coefficients \(f^​\(S\)\\hat\{f\}\(S\)\)\.This heterogeneous support enables\\afsto handle problems with native PB formulations directly, avoiding the representational blowup of CNF translation\.

### 3\.3Gradient Descent and other search Algorithms

We followFastFourierSATand select Projected Gradient Descent \(PGD, Algorithm[1](https://arxiv.org/html/2606.06641#algorithm1)\) as a fast baseline search algorithm for\\afs\. We also provide various other algorithms that possess compatible implementations\. For any selected search algorithm, we take a vector\-map across a batch ofBBcandidate assignments\. This is realised byJAXin GPU warps, whereby every candidate assignment has dedicated memory and a streaming processor, and the chosen algorithm executes in lockstep across every candidate in the batch\. For the choice of PGD, each descent operates numerically independently: a line search determines the step size, a gradient step is taken, and the result is projected back onto the bounded subspace𝒬n\\mathcal\{Q\}^\{n\}if we happen to step out\. Termination occurs when either a convergence criterion is met \(e\.g\. new location is withinε\\varepsilonof the previous location\) or given step\-limit is reached\. The convergence criteria is checked after every step, and some candidates in the batch may converge sooner than others\. As parallelism is achieved via warps, the converged candidates execute the equivalent of no\-ops until the entire batch finishes\.

Input:Fourier expanded formula

FEϕ\\texttt\{FE\}\_\{\\phi\}, initial valuation

X\(0\)\\textbf\{X\}\_\{\(0\)\}, maximum iterations

dd, convergence threshold

δ\\delta, maximum step\-size

ss, bounds

𝒬\\mathcal\{Q\}
Output:Finishing assignment

X¯\\bar\{\\textbf\{X\}\}, evaluation

FEϕ​\(X¯\)\\texttt\{FE\}\_\{\\phi\}\(\\bar\{\\textbf\{X\}\}\), unsatisfied clause count

\#​\(∑m1​\.FEgm​\(X¯\)\>0\)\\\#\\left\(\\sum\_\{m\}1\\text\{ \. \}\\texttt\{FE\}\_\{g\_\{m\}\}\(\\bar\{\\textbf\{X\}\}\)\>0\\right\), iterations

tt
for*t←1t\\leftarrow 1todd*do

𝜼←lineSearch​\(FEϕ,∇\(FEϕ\),s\)\\bm\{\\eta\}\\leftarrow\\text\{lineSearch\}\\left\(\\texttt\{FE\}\_\{\\phi\},\\nabla\(\\texttt\{FE\}\_\{\\phi\}\),s\\right\);

//Determine step length

X\(t\)←X\(t−1\)−𝜼⋅∇\(FEϕ​\(X\(t−1\)\)\)\\textbf\{X\}\_\{\(t\)\}\\leftarrow\\textbf\{X\}\_\{\(t\-1\)\}\-\\bm\{\\eta\}\\cdot\\nabla\\left\(\\texttt\{FE\}\_\{\\phi\}\(\\textbf\{X\}\_\{\(t\-1\)\}\)\\right\);

//Take descent step

X\(t\)←projectToBounds​\(X\(t\),𝒬\)\\textbf\{X\}\_\{\(t\)\}\\leftarrow\\text\{projectToBounds\}\(\\textbf\{X\}\_\{\(t\)\},\\mathcal\{Q\}\);

//Return to bounds

i​s​\_​s​a​t←Trueis\\\_sat\\leftarrow\\texttt\{True\}if

unsat\_count​\(ϕ,X\(t\)\)=0\\text\{unsat\\\_count\}\(\\phi,\\textbf\{X\}\_\{\(t\)\}\)=0elseFalse;

if*η<δ​or​i​s​\_​s​a​t\\eta<\\delta\\text\{ \{or\} \}is\\\_sat*then//Converged or SAT

return

X\(t\),FEϕ​\(X\(t\)\),unsat\_count​\(ϕ,X\(t\)\),t\\textbf\{X\}\_\{\(t\)\},\\texttt\{FE\}\_\{\\phi\}\(\\textbf\{X\}\_\{\(t\)\}\),\\text\{unsat\\\_count\}\(\\phi,\\textbf\{X\}\_\{\(t\)\}\),t;

end if

end for

return

X\(t\),FEϕ​\(X\(t\)\),unsat\_count​\(ϕ,X\(t\)\),d\\textbf\{X\}\_\{\(t\)\},\\texttt\{FE\}\_\{\\phi\}\(\\textbf\{X\}\_\{\(t\)\}\),\\text\{unsat\\\_count\}\(\\phi,\\textbf\{X\}\_\{\(t\)\}\),d;

Algorithm 1Projected Gradient Descent
### 3\.4Partial Assignment Support

\\afs

accepts a partial variable assignment as input, fixing specified variables across all batched starting valuations while randomising the remainder\.JAXsupports gradient masking \(zeroing\), which forces zero gradient on the relevant variables during the automatic differentiation pass resulting in a zero step length\. This feature enables integration into decomposition frameworks where a partial assignment from systematic search is completed by CLS, where some orchestrator may track candidate assignments for solvers in a portfolio to improve or advise upon\.\\afsevenly divides up a given batch of candidate assignment amongst multiple partial assignment if they are provided\. During clause processing,\\afswill also detect trivial unit literals and either fold them into all provided partial assignments, or convert to a universal partial assignment\.

### 3\.5Floating\-Point Limits on Constraint Length

The DFT\-based evaluation for a constraint of lengthnninvolves products∏i=1n\(ωj\+xi\)\\prod\_\{i=1\}^\{n\}\(\\omega^\{j\}\+x\_\{i\}\)which, for rootsωj≈1\\omega^\{j\}\\approx 1and literal valuesxi≈1x\_\{i\}\\approx 1, can reach magnitudes of2n2^\{n\}, while forωj≈−1\\omega^\{j\}\\approx\-1these products approach values on the order ofεn\\varepsilon^\{n\}where\|ε\|→0\\left\\lvert\\varepsilon\\right\\rvert\\to 0\. For IEEE\-754 64\-bit arithmetic with machine epsilonϵ=2−52\\epsilon=2^\{\-52\}, the dynamic range of these terms exceeds representable precision whenn⪆50n\\gtrapprox 50, causing catastrophic cancellation in the inverse DFT, producing incorrect and impermissible evaluations\|FEϕ\|≫1\|\\texttt\{FE\}\_\{\\phi\}\|\\gg 1\.

We empirically observe degenerate solver behaviour—exploding and vanishing gradients—for constraints of lengthn≥48n\\geq 48for which many of variables trend toward the same truth value \(e\.g\. all but one variable assignment in anAMOconstraint of length 50 will tend to false \(1\)\)\. This establishes a practical ceiling for CLS on current GPU architectures that lack practical support for extended\-precision arithmetic\.

#### 3\.5\.1Tailored DFT Construction

Standard DFT implementations \(such as those from common scientific processing packages e\.g\., fromSciPy\[2020SciPy\-NMeth\]\) introduce cumulative errors through repeated exponentiation of roots of unity, breaking the conjugate symmetry required for precise cancellation in the inverse DFT\. Since\\afs’s algebraic foundation is sensitive to these inaccuracies, we implement or own tailored procedures:

- •Forward and inverse DFT matrices are constructed to guarantee exact conjugate symmetry between paired terms\.
- •Combinatorial terms for Fourier coefficients are computed using arbitrary\-precision integer arithmetic, with division and conversion to float deferred to the latest possible stage\.
- •Where algebraic symmetries exist, terms are explicitly mirrored rather than independently computed\.

## 4Performance Evaluation

All experiments were conducted on theGadisupercomputer \(NCI, Australia\), using nodes from the Volta GPU partition\. Each node provides four NVIDIA Tesla V100 GPUs \(32GB HBM2 each,≈\\approx900GB/s bandwidth\), Intel Xeon Cascade Lake CPUs \(48 cores\), and≈\\approx192GB system RAM\.

### 4\.1Comparison withFastFourierSAT

![Refer to caption](https://arxiv.org/html/2606.06641v1/x1.png)
![Refer to caption](https://arxiv.org/html/2606.06641v1/x2.png)

Figure 1:Replicated benchmarks for random cardinality*\(Subfig\. a\)*and parity learning \(xor\) with errors*\(Subfig\. b\)*as perFastFourierSAT\[cen2025massively\], with PAR\-2 scores and timeouts indicated\. We run\\afsconfigured as close to possible toFastFourierSAT, replicating the batch size choice as indicated in the original benchmark\. We also run\\afsin*max\-through*mode, where the batch size is selected to target maximum efficiency \(see §[4\.2](https://arxiv.org/html/2606.06641#S4.SS2)\)\.We replicate the cardinality constraint and parity\-learning benchmarks from Cen et al\.\[cen2025massively\]\. Figure[1](https://arxiv.org/html/2606.06641#S4.F1)compares cumulative solution times\. We apply certainJAXoptimisations set for\\afstoFastFourierSATto provide a more level playing field, and note that the parity learning results shown forFastFourierSATare better than those presented in their paper due to bugs in the published code we have corrected\. Key observations:

- •\\afs achieves equivalent performance in the worst\-case, and otherwise universally improves uponFastFourierSAT, demonstrating that our optimisations strictly improved upon the proof\-of\-concept on comparable benchmarks\.
- •\\afs exhibits an obviously lower baseline times thanFastFourierSAT, indicating more efficient pre\-processing, GPU kernel optimisation/compilation, and general overheads\. For fairness, we take the best time of several runs for each to allow for disk latency and compilation caches to be populated\.
- •\\afs consumes substantially less GPU memory by storing only minimal auxiliary data \(e\.g\. highly optimised DFT matrices, clause and literal arrays, etc\.\) and computing closures over this data, enabling the compilers to better optimise which subsequently enables larger batch sizes\.

### 4\.2GPU Memory and Throughput Characteristics

During testing it was noticed that as total GPU memory utilisation increased \(reflecting larger batch sizes\), the total number of completed gradient descents decreased both relatively \(with respect to the amount of memory consumed\) and in some cases absolutely\. We investigate the total number of completed searches we can take per unit time, a measure we call*throughput*, as a function of GPU memory consumption and GPUs available across several problem domains \(Figure[2](https://arxiv.org/html/2606.06641#S4.F2)\)\. We observe clearly that throughput peaks when the memory used is approximately0\.1%0\.1\\%–1%1\\%of total GPU memory consumed\. When exploring this effect across various GPUs, we find this percentage varies, but upon closer inspection the value is more tightly correlated twice the cache memory available on the device\. We observe a logarithmic decay \(Figure[2](https://arxiv.org/html/2606.06641#S4.F2)b\) in throughput beyond this peak, rather than the monotonic increases one might expect when parallelising a problem up to memory saturation\.

![Refer to caption](https://arxiv.org/html/2606.06641v1/x3.png)
![Refer to caption](https://arxiv.org/html/2606.06641v1/x4.png)

Figure 2:Throughput vs measured memory consumption across various problem classes\.
*\(Subfig\. a\)*shows near linear scaling when parallelising across several accelerators on various hard random 3SAT problems, with the occasional exception due to properties of individual test cases\.
*\(Subfig\. b\)*shows scaled throughput for various benchmarking problems, confirming the dynamics of throughput vs\. batch size \(memory consumption\)\. A polynomial trendline peaks between0\.1%0\.1\\%and1%1\\%\.Across various domains we also observe the optimal batch size \(and subsequently peak throughput\) is governed primarily by maximum constraint length rather than the absolute number of variables or constraints\. This is dominated ultimately by the requirement to computeO​\(k2\)O\(k^\{2\}\)DFT matrices forkklength clauses\. As the size of the DFT grows quadratically it rapidly cuts the amount of memory usable for the remainder of the algorithm, resulting in lower batch sizes\. We hypothesise that the throughput peak corresponds to saturation of low\-level processing core caches rather than main memory, whereby optimal pipelining is achieved with next to no cache invalidation or page\-faulting\. Supporting this, the ratio of total cache to total memory on the V100 architecture is on the order of0\.1%0\.1\\%—consistent with the observed peak locations\.

This finding has practical implications: for problems with long constraints, smaller batch sizes yield higher throughput than larger ones that exceed cache capacity\. For deployment,\\afsprovides batch\-size tuning that balances throughput against total search coverage\. As indicated in §[3\.5](https://arxiv.org/html/2606.06641#S3.SS5), long constraints also degrade numerical stability\. The combination of both of these effects suggests that future improvements will look to trade\-off the compactness of single constraints for the efficiency gains of constraint decompositions\.

JAX’s sharding mechanism distributes batched valuations across available GPUs in a compute\-follows\-data paradigm\. Since gradient descents are independent across each assignment in the batch, communication overhead is minimal\. We observe near\-linear scaling in throughput with increasing GPU count across all tested problem domains and Figure[2](https://arxiv.org/html/2606.06641#S4.F2)a demonstrates this for hard random 3SAT problems\. The absence of an impacting overhead cost confirms the suitability of CLS for multi\-accelerator deployment\.

## 5Conclusion and Future Work

\\afs

advances the state of CLS\-based SAT solving from proof\-of\-concept to a practical, extensible tool\. By engineering support for heterogeneous pseudo\-Boolean constraints, achieving substantial performance improvements over the baseline, and demonstrating scalable multi\-GPU execution,\\afsestablishes a foundation for GPU\-accelerated SAT solving in accelerator\-rich environments\. The identification of precision\-driven constraint\-length limits and cache\-driven throughput characteristics provides actionable guidance for deployment\. Future development will target adaptive constraint weighting, integration with systematic solvers via decomposition frameworks, and exploitation of emerging higher\-precision accelerator arithmetic\.

##### Constraint length ceiling and memory consumption\.

As established in §[3\.5](https://arxiv.org/html/2606.06641#S3.SS5)and §[4\.2](https://arxiv.org/html/2606.06641#S4.SS2), floating\-point precision limits constraint length ton≈50n\\approx 50on 64\-bit hardware and dominates memory consumption\. Problems with naturally longer constraints \(e\.g\., global cardinality constraints in large graph colouring\) can be addressed by decomposing these into shorter equivalent constraints\. For all but EK and CARD, there are simple linear decompositions toO​\(k/50\)O\(k/50\)equivalent PB constraints of the same type, preserving some of the benefits of native PB representation\.

##### Incompleteness\.

\\afs

is an incomplete solver: it cannot prove unsatisfiability, as doing so would require solving the global optimisation problem to certifiable optimality\. However,\\afsfunctions naturally as a MaxSAT solver, providing best\-effort solutions with unsatisfied\-constraint counts\. Unbounded global optimisation methods are a potential consideration, and methods involving Moreau envelopes\[cen2025massively\]or unconstrained formulations with extreme penalties\[zhang2025thinkingboxhybridsat\]\.

##### Second\-order methods\.

The multi\-linear objective is thoroughly populated with saddle points, at which first\-order PGD can stall\. Bounded second\-order methods \(e\.g\., L\-BFGS\-B\) would help escape saddle points, but existing implementations in theJAXecosystem—notablyJAXOpt\[jaxopt\_implicit\_diff\]andOptimistix\[optimistix2024\]—do not correctly respect bounds or are not production\-ready\. Higher order methods will necessarily decrease throughput due to memory required to compute and store higher order gradients\.

## References

## Appendix AInput Format Specification

\\afs

accepts problem instances in standard DIMACS CNF or an extended DIMACS\-like hybrid format supporting heterogeneous symmetric pseudo\-Boolean constraints\.

##### DIMACS\-like PB constraint grammar

The hybrid PB problem grammar is defined as follows, where tokens are whitespace\-separated:

`The type identifiers map to constraint semantics as follows:`

`Identifiers Type Semantics \(none\) CNF Standard disjunction x, xor XOR Odd number of literals true n, nae NAE Not all literals equal a, amo AMO At most one literal true e, eo EO Exactly one literal true k, ek EK Exactly kk literals true d, card CARD Cardinality threshold Table 2: Constraint type identifiers and semantics\.``For CARD constraints with a plain integer threshold \(no operator\): positive kk defaults to ≥k\\geq k; negative kk defaults to <\|k\|<\|k\|\. All CARD constraints are internally normalised to ≥\\geq or << forms\. Examples A\.1 Automatic Simplification The parser applies the following reductions: • Unit EO and CNF constraints are extracted as unit\-propagated prefix assignments\. • Unit AMO constraints are discarded \(trivially satisfied\)\. • CARD\-11 constraints are reduced to CNF; EK\-11 to EO\. • CARD\-nn and EK\-nn \(where kk equals the number of literals\) are converted to unit prefix assignments\. • CARD\-0 constraints are discarded \(trivially satisfied\); EK\-0 constraints yield negated prefix assignments\. • Conflicting unit literals, if detected, immediately raise an unsatisfiability error\.`