\section{Introduction}
Constraint satisfaction problems (CSPs) are widely used in AI, with applications in optimization, control, and planning \cite{russell2003artificial}. 
While many classes of CSPs are intractable in the worst case~\cite{cook1971complexity}, many real-world CSP instances are easy to solve in practice~\cite{vardi2014boolean}.  
As a result, there has been significant interest in understanding the average-case complexity of CSPs from multiple communities, such as AI, theoretical computer science, physics, and combinatorics~\cite{biere2009handbook}. 

Random CSPs are an important model for studying the average-case complexity of CSPs. Past works have proposed several distributional models for random CSPs \cite{molloy2003models,creignou2003generalized}. An interesting feature arising from many models is a phase transition phenomenon that occurs as one changes the ratio of number of constraints, $m$, to the number of variables, $n$. Empirical results \cite{mitchell1992hard} show that for many classes of CSPs, randomly generated instances are satisfiable with probability near 1 when $m/n$ is below a certain threshold. For $m/n$ larger than a threshold, the probability of satisfiability is close to 0. The statistical physics, computer science, and mathematics communities have focused much attention on identifying these threshold locations \cite{achlioptas2005rigorous,biere2009handbook}. 

Phase transitions for common classes of CSPs such as $k$-SAT and $k$-XORSAT are very well-studied. For $k$-SAT, researchers struggled to find tight lower bounds on the satisfiability threshold until the breakthrough work of \citeauthor{achlioptas2002asymptotic}, which provided lower bounds a constant factor away from the upper bounds. Later works closed this gap for $k$-SAT \cite{coja2013going,coja2016asymptotic}.  More recently, \citeauthor{Dudek2016CombiningTK} \shortcite{Dudek2016CombiningTK} also studied the satisfiability threshold for a more general CSP class, namely $k$-CNF-XOR, where both $k$-SAT and $k$-XORSAT constraints can be used. The results and analyses from these works, however, are all specific to the constraint classes studied.

In this paper, we provide new lower bounds on the location of the satisfiability threshold that hold for general boolean CSP classes. We focus on the setting where CSPs are generated by a single constraint type, though our analysis can extend to the setting with uniform mixtures of different constraint functions. We extend techniques from \cite{achlioptas2004threshold} and build on \cite{creignou2003generalized}, which proposes a distributional model for generating CSPs and provides lower bounds on the satisfiability threshold for these models. The significance of our work is that \textit{our bounds hold for all functions that could be used to generate random CSP instances}. The lower bounds from \cite{creignou2003generalized} are also broadly applicable, but they are looser than ours because they do not depend on constraint-specific properties. Our lower bounds are often tight (within a constant factor of upper bounds for many CSP classes) because they depend on specific properties of the Fourier spectrum of the function used to generate the random CSPs. Since these properties are simple to compute for any constraint function, our lower bounds are broadly applicable too. 

The Fourier analysis of boolean functions \cite{o2014analysis} will be vital for obtaining our main results. Expressing functions in the Fourier basis allows for clean analyses of random constraints \cite{friedgut1999sharp,barak2015beating,achim2016beyond}. Our use of Fourier analysis is inspired by the work of \citeauthor{achim2016beyond} \shortcite{achim2016beyond}, who analyze the Fourier spectra of random hash functions used as constraints in CSP-based model counting. We show that the Fourier spectrum of our constraint-generating function controls the level of spatial correlation in the set of satisfying assignments to the random CSP. If the Fourier spectrum is concentrated on first and second order coefficients (corresponding to ``low frequencies''), this correlation will be very high, roughly increasing the variance of the number of solutions to a random CSP and decreasing the probability of satisfiability. In related work, \citeauthor{montanari2011reconstruction} \shortcite{montanari2011reconstruction} also use Fourier analysis to provide tight thresholds in the case where odd Fourier coefficients are all zero. 

\section{Notation and Preliminaries}
\label{sec:notation}
In this section, we will introduce the preliminaries necessary for presenting our main theorem. First, we formally define our distribution for generating random CSPs, inspired from \cite{creignou2003generalized}. 

We will use $n$ and $m$ to denote the number of variables and number of constraints in our CSPs, respectively. We will also let $f : \{-1, 1\}^k \rightarrow \{0, 1\}$ denote a binary function and refer to $f$ as our constraint function. Often we use the term ``solution set of $f$'' to refer to the set $\{u : f(u) = 1\}$. Using the constraint function $f$, we create constraints by applying $f$ to a signed subset of $k$ variables. 

\begin{definition} [Constraint]
Let $I = (i_1, \ldots, i_k)$ be an ordered tuple of $k$ indices in $[n]$, and let $s$ be a sign vector from $\{-1, 1\}^k$. Given a vector $\sigma \in \{-1, 1\}^n$, we will define the vector $\sigma_{I, s}$ of size $k$ as follows: 
\begin{equation}
\label{eq:sigmaIs}
\sigma_{I, s} = (s_1\sigma_{i_1}, \ldots, s_k\sigma_{i_k})
\end{equation}
Now we can denote the application of $f$ to these indices by $f_{I, s}(\sigma) = f(\sigma_{I, s})$. We call $f_{I, s}$ a constraint, and we say that $\sigma \in \{-1, 1\}^n$ satisfies the constraint $f_{I, s}$ if $f_{I, s}(\sigma) = 1$.
\end{definition}
The definition of a CSP generated from $f$ follows. 
\begin{definition} [CSP generated from $f$]
We will represent a CSP with $m$ constraints and $n$ variables generated from $f$ as a collection of constraints $C_f(n, m) = \{f_{I_1, s_1}, \ldots, f_{I_m, s_m}\}$. Then $\sigma \in \{-1, 1\}^n$ satisfies $C_f(n, m)$ if $\sigma$ satisfies $f_{I_j, s_j}$ for $j = 1, \ldots, m$.
\end{definition}
\begin{example} [3-SAT]
Let $f : \{-1, 1\}^3 \rightarrow \{1, 0\}$ where $f(u) = 0$ for $u = (-1, -1, -1)$ and $f(u) = 1$ for all other $u$. Then $f$ is the constraint function for the 3-SAT problem.
\end{example}
With these basic definitions in place, we are ready to introduce the model for random CSPs. 
\subsection{Random CSPs}
\label{sec:randcsp}
We discuss our model for randomly generating  $C_f(n, m)$, and formally define a ``satisfiability threshold.''

To generate instances of $C_f(n, m)$, we simply choose $I_1, \ldots, I_m$ and $s_1, \ldots, s_m$ uniformly at random. For completeness, we sample without replacement, i.e. there are no repeated variables in a constraint, and no duplicate constraints. However, sampling with replacement does not affect final results. For the rest of this paper we abuse notation and let $C_f(n, m)$ denote a randomly generated CSP instance following this model. 

Now we formally discuss satisfiability thresholds. We let $r = m/n$. For many constraint functions $f$, there exist thresholds $r_{f,\text{sat}}$ and $r_{f,\text{unsat}}$ such that 
\begin{align*}
\lim_{n \rightarrow \infty} \Pr[C_f(n, rn) \ \text{is satisfiable}] = 
\begin{cases}
1 \ \text{if} &\ r < r_{f,\text{sat}}\\
0 \ \text{if} &\ r > r_{f,\text{unsat}}
\end{cases}
\end{align*}
In general, it is unknown whether $r_{f,\text{sat}} = r_{f,\text{unsat}}$, but for some problems such as $k$-SAT (for large $k$) and $k$-XORSAT, affirmative results exist \cite{ding2015proof,pittel2016satisfiability}. If $r_{f,\text{sat}} = r_{f,\text{unsat}}$, then we say that the random CSP $C_f(n, m)$ exhibits a sharp threshold in $m/n$. 

We are concerned with finding lower bounds $r_{f,\text{low}}$ on $r_{f,\text{unsat}}$ such that there exists a constant $C > 0$ independent of $n$ so that for sufficiently large $n$, 
\begin{equation}
\label{eq:rflower}
\Pr[C_f(n, rn) \ \text{is satisfiable}] > C \ \text{for all} \ r < r_{f,\text{low}}
\end{equation}
For CSP classes with a sharp threshold, \eqref{eq:rflower} implies that 
\begin{align*}
\lim_{n \rightarrow \infty} \Pr[C_f(n, rn) \ \text{is satisfiable}] = 1 \ \text{for all} \ r < r_{f,\text{low}}
\end{align*}
We also wish to find upper bounds $r_{f,\text{up}}$ such that 
\begin{align*}
\lim_{n \rightarrow \infty} \Pr[C_f(n, rn) \ \text{is satisfiable}] = 0 \ \text{for all} \ r > r_{f,\text{up}}
\end{align*}
For an example of these quantities instantiated on a concrete example, refer to the experiments in Section \ref{sec:experiments}.

We provide a value for $r_{f,\text{up}}$ which was derived earlier in \cite{creignou2003generalized}. \citeauthor{dubois2001upper} \shortcite{dubois2001upper} and \citeauthor{creignou2007expected} \shortcite{creignou2007expected} provide methods for obtaining tighter upper bounds, but the looser values that we use are sufficient for showing that $r_{f,\text{low}}$ is on the same asymptotic order as $r_{f,\text{unsat}}$ for many choices of $f$. 

The bound $r_{f,\text{low}}$ depends on both the symmetry and size of the solution set of $f$. The more assignments $u \in \{-1, 1\}^k$ such that $f(u) = 1$, the more likely it is that each constraint is satisfied. Increased symmetry reduces the variance in the number of solutions to $C_{f}(n, rn)$, so solutions are more spread out among possible CSPs in our class and the probability that $C_f(n, rn)$ will have a solution is higher. We will formally quantify this symmetry in terms of the Fourier spectrum of $f$, which we introduce next. 

\subsection{Fourier Expansion of Boolean Functions}
We discuss basics of Fourier analysis of boolean functions. For a detailed review, refer to \cite{o2014analysis}. We define the vector space $\mathcal{F}_k$ of all functions mapping $\{-1, 1\}^k$ to $\mathbb{R}$. The set $\mathcal{F}_k$ has the inner product $\langle f_1, f_2 \rangle = \sum_{u \in \{-1, 1\}^k} f_1(u)f_2(u)/2^k$ for any $f_1, f_2 \in \mathcal{F}_k$.

This inner product space has orthonormal basis vectors $\chi_S$, where the parity functions $\chi_S$ follow $\chi_S(u) = \prod_{i \in S} u_i$ for all $S \subseteq [k]$, subsets of the $k$ indices. Because $(\chi_S)_{S \subseteq [k]}$ forms an orthonormal basis, if we write 
\begin{equation}
\label{eq:hatfS}
\hat{f}(S) = \langle f, \chi_S \rangle = \frac{1}{2^k} \sum_{u \in \{-1, 1\}^k} f(u)\chi_S(u)
\end{equation}
then we can write $f$ as a linear combination of these vectors: $f = \sum_{S \subseteq [k]} \hat{f}(S) \chi_S$. We note that when $S = \emptyset$, the empty set, $\hat{f}(\emptyset)$ is simply the average of $f$ over $\{-1, 1\}^k$. We will refer to the coefficients $(\hat{f}(S))_{S \subseteq [k]}$ as the Fourier spectrum of $f$. Since these coefficients are well-studied in theoretical computer science \cite{o2014analysis}, the Fourier coefficients of many boolean functions are easily obtained.

\begin{example} [3-SAT] 
\label{ex:3satfourier}
For $3$-SAT, $\hat{f}(\emptyset) = 7/8$. $\hat{f}(\{1\}) = \hat{f}(\{2\}) = \hat{f}(\{3\}) = 1/8$, and $\hat{f}(\{1, 2\}) = \hat{f}(\{2, 3\}) = \hat{f}(\{1, 3\}) = -1/8$.
\end{example}

Representing $f$ in the Fourier bases will facilitate our proofs, providing a simple way to express expectations over our random CSPs. The Fourier spectrum can also provide a measure of ``symmetry'' in $f$ - if some values of $\hat{f}(S)$ are high where $|S| = 1$, then satisfying assignments to $f$ are more skewed in the variable corresponding to $S$. We will show how this impacts satisfiability in Section \ref{sec:boundexplained}.

\begin{figure*}
\centering
\begin{tabular}{|c | c | c| c|}
\hline
\Tstrut
CSP class ($f$) & Best lower bound on $r_{f,\text{unsat}}$ & Our bound $r_{f,\text{low}}$ & Upper bound $r_{f,\text{up}}$\\ \hline
\Tstrut
$k$-XORSAT & 1 & $\frac{1}{2}$ & 1\\
\Tstrut
$k$-SAT & $2^k\ln 2 - \frac{1 + \ln 2}{2} - o_k(1)$ & $2^{k - 1} - O(k)$ & $2^k \ln 2$\\ 
\Tstrut
$k$-NAESAT & $2^{k - 1} \ln 2 - \frac{\ln 2}{2} - \frac{1}{4} - o_k(1)$ & $2^{k - 2} - \frac{1}{2}$ & $2^{k - 1}\ln 2$ \\
\Tstrut
$k$-MAJORITY & ? & $\frac{\frac{1}{2} - k\binom{k - 1} {\frac{k - 1}{2}}^2 2^{-2k +1}}{1 + k\binom{k - 1}{\frac{k - 1}{2}}^2 2^{-2k+ 2}} = 0.111 - o_a(1)$ & $1$  \\
\Tstrut
$a$-MAJ $\otimes 3$-MAJ & ? & $\frac{\frac{1}{2} - 3a\binom{a - 1}{\frac{a - 1}{2}}^2 2^{-2a - 1}}{1 + 3a\binom{a - 1}{\frac{a - 1}{2}}^2 2^{-2a - 2}} = 0.177 - o_a(1)$ & $1$\\
\Tstrut
$k$-MOD-3 & ? & $\frac{1}{4} - o_k(1)$ & $\frac{\ln 2}{\ln 3} + o_k(1)$\\
\Tstrut
OR$_b \otimes$ XOR$_a$ & ? & $2^{b - 1} - 1/2$ & $2^{b - 1} \ln 2$\\
\hline
\end{tabular}
\caption{We compare the best known lower bounds on the satisfiability threshold to our lower and upper bounds. For $k$-XORSAT \cite{pittel2016satisfiability}, $k$-SAT \cite{ding2015proof}, and $k$-NAESAT \cite{coja2012catching}, the numbers listed are known as exact sharp threshold locations. For the last four, we do not know of existing lower bounds. $\otimes$ is the composition operator for boolean functions, and we define these functions in Section \ref{sec:figfunc}.}
\label{fig:bounds}
\end{figure*}

\section{Main Results}
\label{sec:mainresult}
We provide simple formula for $r_{f,\text{up}}$. A similar result is in \cite{creignou2003generalized}, and the full proof is in the appendix.
\begin{proposition}
\label{prop:upper}
For all constraint functions $f$, let 
\begin{align*}
r_{f,\text{up}} = \frac{\log 2}{\log 1/\hat{f}(\emptyset)} < \frac{\log 2}{1 - \hat{f}(\emptyset)}
\end{align*}
If $r \ge r_{f,\text{up}}$, $\lim_{n \rightarrow \infty} \Pr[C_f(n, rn) \ \text{is satisfiable}] = 0$. 
\end{proposition}
\begin{proof} [Proof Sketch]
We compute the expected solution count for $C_f(n, rn)$. The expected solution count will scale with $\hat{f}(\emptyset)$, since  $2^k \hat{f}(\emptyset)$ is simply the number of $u \in \{-1, 1\}^k$ where $f(u) = 1$ and therefore governs how easily each constraint will be satisfied. If $r > r_{f,\text{up}}$, the expected solution count converges to 0 as $n \rightarrow \infty$, so Markov's inequality implies that the probability that a solution exists goes to 0.
\end{proof}

Next, we will present our value for $r_{f,\text{low}}$. First, some notation: let $U = \{u \in \{-1, 1\}^k : f(u) = 1\}$, and let $A$ be the $k \times |U|$ matrix whose columns are the elements of $u$. We will use $A^+$ to denote the Moore-Penrose pseudoinverse of $A$. For a reference on this, see \cite{barata2012moore}. Finally, let $\boldsymbol{1}$ be the $|U|$-dimensional vector of 1's.
\begin{example} [$3$-SAT] For $3$-SAT, $k = 3$ and $|U| = 7$, and we can write $A$ as follows (up to permutation of its columns):
\begin{align*}
A = 
\left[
\begin{array} {r r r r r r r}
1 & 1 & 1 & 1 & -1 & -1 & -1\\
1 & 1 & -1 & -1 & 1 & 1 & -1\\
1 & -1 & 1 & -1 & 1 & -1 & 1
\end{array}
\right]
\end{align*}
where columns of $A$ satisfy the $3$-SAT constraint function.
\end{example}
The following main theorem provides the \textit{first computable equation for obtaining lower bounds that are specific to the constraint-generating function $f$.} 

\begin{theorem}
\label{thm:main}
For all constraint functions $f$, let 
\begin{align*}
r_{f,\text{low}} = \frac{1}{2} \frac{c}{1 - c} \text{ where } c = \hat{f}(\emptyset) - \frac{\boldsymbol{1}^T A^+ A \boldsymbol{1}}{2^k}
\end{align*}
If $r < r_{f,\text{low}}$, then there exists a constant $C>0$ such that $\lim_{n \rightarrow \infty} \Pr[C_f(n, rn) \ \text{is satisfiable}] > C$. 
\end{theorem}

The lower bound $r_{f,\text{low}}$ is an increasing function of $c$ which is dependent on two quantities. First, with higher values of $\hat{f}(\emptyset)$, $C_f(n, rn)$ will have more satisfying assignments on average, so $c$ and the threshold value will be higher. Second, $c$ depends on the level of symmetry in the solution set of $f$, which we will show is connected to the Fourier spectrum of $f$. We explain this dependence in Section \ref{sec:boundexplained}. 

In comparison, \citeauthor{creignou2003generalized} \shortcite{creignou2003generalized} obtain lower bounds which depend only on the arity of $f$. While we cannot make an exact comparison because \citeauthor{creignou2003generalized} use a different random CSP ensemble, for reference, they provide the general lower bound of $1/(ke^k - k)$ expected constraints per variable for functions of arity $k$. Our bounds are much tighter because of their specificity while remaining simple to compute. To demonstrate, we instantiate our bounds for some example constraint functions in Figure \ref{fig:bounds}. Whereas their bounds are exponentially decreasing in $k$, our bounds are constant or increasing in $k$ for the functions shown. 

\subsection{Constraint Functions in Figure \ref{fig:bounds}}
\label{sec:figfunc}
We define the constraint functions in Figure \ref{fig:bounds}. Unless specified otherwise, they will be in the form $f : \{-1, 1\}^{k} \rightarrow \{0, 1\}$. 
\begin{enumerate}
\item $k$-SAT: $f(u) = 0$ if $u$ is the all negative ones vector, and $f(u) = 1$ otherwise.
\item $k$-XORSAT: $f(u) = \mathbbm{1}(\chi_{[k]}(u) = -1)$
\item $k$-NAESAT: $f(u) = 0$ if $u$ is the all negative ones or all ones vector, and $f(u) = 1$ otherwise.
\item $k$-MAJORITY: Defined when $k$ is odd, $f(u) = 1$ if more than half of the variables of $u$ are $1$.
\item $a$-MAJ $\otimes$ $3$-MAJ: Defined when $a$ is odd, where $f : \{-1, 1\}^{3a} \rightarrow \{0, 1\}$. Defined as the composition of $a$-MAJORITY on $a$ groups of $3$-MAJORITY, as follows: 
\begin{align*}
f(u_1, \ldots, u_{3a}) =\\ f_{a-\text{MAJ}}(f_{3-\text{MAJ}}(u_1, u_2, u_3), \ldots, \\ f_{3-\text{MAJ}}(u_{3a - 2}, u_{3a -1}, u_{3a}))
\end{align*}
\item $k$-MOD-3: $f(u) = 1$ when the number of 1's in $u$ is divisible by $3$, and 0 otherwise.
\item OR$_b \otimes$ XOR$_a$: In this case, $f : \{-1, 1\}^{ab} \rightarrow \{0, 1\}$, and $f$ is the composition of a OR over $b$ groups of XORs over $a$ variables, as follows: 
\begin{align*}
f(u_{1}, \ldots, u_{ab}) =\\
f_{\text{OR}_b}(f_{\text{XOR}_a}(u_1, \ldots, u_a), \ldots,\\
f_{\text{XOR}_{a}}(u_{ab - a + 1}, \ldots, u_{ab}))
\end{align*}
\end{enumerate}

While the last four constraint functions have not been analyzed much in the existing CSP literature, these types of general constraints are of practical interest because of \cite{achim2016beyond}, which performs probabilistic inference by solving CSPs based on arbitrary hash functions. For example, \citeauthor{achim2016beyond} \shortcite{achim2016beyond} show that MAJORITY constraints are effective in practice for solving probabilistic inference problems. 
\subsection{Connecting Bounds with Fourier Spectrum}
\label{sec:boundexplained}
We explain how the Fourier spectrum can help us interpret Theorem \ref{thm:main}. We first show the connection between $c$ and the Fourier spectrum. 
Let $\hat{f}_{S : |S| = 1}$ be the $k$-dimensional vector whose entries are Fourier coefficients of $f$ for size 1 sets. Let $B$ be the $k \times k$ matrix with diagonal entries $B_{ii} = \hat{f}(\emptyset)$ and off-diagonal entries $B_{ij} = \hat{f}(\{i, j\})$ for $i \ne j$. 
\begin{example} [3-SAT]
Following the coefficients in Example \ref{ex:3satfourier}, for $3$-SAT, $\hat{f}_{S : |S| = 1} = (1/8, 1/8, 1/8)$ and 
\begin{align*}
B = \left[
\begin{array} {r r r}
7/8 & -1/8 & -1/8\\
-1/8 & 7/8 & -1/8\\
-1/8 & -1/8 & 7/8
\end{array}
\right]
\end{align*}
\end{example}
\begin{lemma}
\label{lem:fourierc}
When the rows of $A$ are linearly independent, $c = \hat{f}(\emptyset) - \hat{f}_{S : |S| = 1}^T B^{-1} \hat{f}_{S : |S| = 1}$.
\end{lemma}

From this lemma, we see that larger values of $c$ correspond to smaller $\hat{f}_{S : |S| = 1}$. These terms will measure the amount of ``symmetry'' in the solution set for $f$. 
The matrix $B$ and vectors $\hat{f}_{S : |S| = 1}$ are easily obtained for many $f$ since Fourier coefficients are well-studied \cite{o2014analysis}.  

Figure \ref{fig:bounds} shows how $\hat{f}(\emptyset) - c$ and $f$ relate. Since $k$-SAT has a mostly symmetric solution set, $\hat{f}_{k-\text{SAT}}(\emptyset) - c = O(k/2^{2k})$ since $c = 1 - 2^{-k} - O(k/2^{2k})$, which is small compared to $\hat{f}_{k-\text{SAT}}(\emptyset)$. The solution set of $k$-NAESAT is completely symmetric as if $f_{k-\text{NAESAT}}(x) = 1$, then $f_{k-\text{NAESAT}}(-x) = 1$. Thus, $f_{k-\text{NAESAT}}$ has 0 weight on Fourier coefficients for sets with odd size so we can compute that $\hat{f}_{k-\text{NAESAT}}(\emptyset) - c = 0$. $k$-MAJORITY, however, is less symmetric, as shown by larger first order coefficients. Here, the bound in Figure \ref{fig:bounds} gives 
$\lim_{k \rightarrow \infty} \hat{f}_{k-\text{MAJORITY}}(\emptyset) - c = 1/\pi$, which is large compared to $\hat{f}_{k-\text{MAJORITY}}(\emptyset) \approx 1/2$.

\section{Proof Strategy}
Our proof relies on the second moment method, which has been applied with great success to achieve lower bounds for problems such as $k$-SAT \cite{achlioptas2004threshold} and $k$-XORSAT \cite{dubois20023}. The second moment method is based on the following lemma, which can be derived using the Cauchy-Schwarz inequality: 
\begin{lemma}
\label{lem:secondmomentmethod}
Let $X$ be any real-valued random variable. Then
\begin {equation}
\label{eq:secondmom}
Pr[X \ne 0] = \Pr[|X| \ne 0] \ge \frac{\text{E}[|X|]^2}{\text{E}[X^2]} \ge \frac{\text{E}[X]^2}{\text{E}[X^2]}
\end {equation}
\end{lemma}
If $X$ is only nonzero when $C_f(n, rn)$ has a solution, we obtain lower bounds on the probability that a solution exists by upper bounding $\text{E}[X^2]$. For example, we could let $X$ be the number of solutions to $C_f(n, rn)$. However, as shown in \cite{achlioptas2004threshold}, this choice of $X$ fails in most cases. Whether two different assignments satisfy $C_f(n, rn)$ is correlated: if the assignments are close in Hamming distance and one assignment is satisfying, it is more likely that the other is satisfying as well. This will make $\text{E}[X^2]$ much larger than $\text{E}[X]^2$, so \eqref{eq:secondmom} will not provide useful information. Figure \ref{fig:failure} demonstrates this failure for $k$-SAT. \citeauthor{achlioptas2004threshold} \shortcite{achlioptas2004threshold} show formally that the ratio $\text{E}[X]^2/\text{E}[X^2]$ will decrease exponentially (albeit at a slow rate). On the other hand, $k$-NAESAT is ``symmetric'', so the second moment method works directly here. In the plot, $\text{E}[X]^2/\text{E}[X^2]$ for $3$-NAESAT stays above a constant. This also follows formally from our main theorem as well as \cite{achlioptas2002asymptotic}. We formally define our requirements on symmetry in \eqref{eq:derivhalf0}.  

We circumvent this issue by weighting solutions to reduce correlations before applying the second moment method. As in \cite{achlioptas2004threshold}, we use a weighting which factors over constraints in $C_f(n, rn)$ and apply the second moment method to the random variable 
\begin{equation}
\label{eq:X}
X = \sum_{\sigma \in \{-1, 1\}^n} \prod_{c \in C_f(n, rn)}w(\sigma, c)
\end{equation}
where $C_f(n, rn)$ is a collection of constraints $\{f_{I_1, s_1}, \ldots, f_{I_m, s_m}\}$ and the randomness in $X$ comes over the choices of $I_j, s_j$. Now we can restrict our attention to constraint weightings of the form $w(\sigma, f_{I, s}) = w(\sigma_{I, s})$. In the special case where $w(\sigma_{I, s}) = f(\sigma_{I, s})$, $X$ will simply represent the number of solutions to $C_f(n, rn)$. In general, we require $w(\sigma_{I, s}) = 0$ whenever $f(\sigma_{I, s}) = 0$. This way, if $X \ne 0$, then $C_f(n, rn)$ must have a solution.  

For convenience, we assume that the index sets $I_1, \ldots, I_m$ are sampled with replacement. They are chosen uniformly from $[n]^k$. We also allow constraints to be identical. In the appendix, we justify why proofs in this setting carry over to the without-replacement setting in Section \ref{sec:randcsp} and also provide full proofs to the lemmas presented below.  

\begin{figure*}
\centering
\subfloat[$\text{E}{[X]}^2/\text{E}{[X^2]}$ vs. $n$ for $3$-SAT and $3$-NAESAT. $X$ is the solution count, $r = 1$.]
{\includegraphics[width=0.32\textwidth]{naesat_sat_plot_n_max_600_k_3_r_1_num_n_30_n_min_1.png}
\label{fig:failure}
}
\subfloat[$g_w(\alpha)$ when $f$ is $k$-XORSAT and $w = f$.]{\includegraphics[width=0.32\textwidth]{g_w_xor_plot_k_max_12_k_min_4.png}
\label{fig:g_w}}
\subfloat[$\psi_r(\alpha)$ with $f$ as $5$-NAESAT and $w = f$.]{\includegraphics[width=0.32\textwidth]{psi_r_plot_num_r_3_k_5_r_max_13_0_r_min_1_0.png}
\label{fig:psi}}
\caption{For concreteness, we provide sample plots of the relevant quantities in our proofs.}
\end{figure*}

In this setting, we will compute the first and second moments of the $X$ chosen in \eqref{eq:X} in terms of the Fourier spectrum of $w$.
\begin{lemma}
\label{lem:firstmomentsq}
The squared first moment of $X$ is given by 
\begin{equation}
\label{eq:E[X]sq}
\text{E}[X]^2 = 2^{2n}(\hat{w}(\emptyset)^2)^{rn}
\end{equation}
\end{lemma}
\begin{proof}
We can expand $\text{E}[X]$ as follows: 
\begin{align}
\text{E}[X] &= \sum_{\sigma \in \{-1, 1\}^n} \text{E}\left[\prod_{j = 1}^{rn}w(\sigma_{I_j, s_j})\right] \notag \\
&= \sum_{\sigma \in \{-1, 1\}^n} \text{E}[w(\sigma_{I, s})]^{rn} \label{eq:EXint}
\end{align}
where we used the fact that constraints are chosen independently. Now we claim that for any $u \in \{-1, 1\}^k$, 
\begin{align*}
\Pr[\sigma_{I, s} = u] = \frac{1}{2^k}
\end{align*}
This follows from the fact that we choose $s$ uniformly over $\{-1, 1\}^k$ and our definition of $\sigma_{I, s}$ in \eqref{eq:sigmaIs}. Thus, 
\begin{align*}
\text{E}[w(\sigma_{I, s})] &= \sum_{u \in \{-1, 1\}^k} w(u)\Pr[\sigma_{I, s} = u]\\
&= \frac{1}{2^k}\sum_{u \in \{-1, 1\}^k} w(u)\\
&= \hat{w}(\emptyset)
\end{align*}
Plugging back into \eqref{eq:EXint} gives the desired result. 
\end{proof}
Next, we will compute the second moment $\text{E}[X^2]$.
\begin{lemma}
\label{lem:secmoment}
Let $g_w(\alpha) = \sum_{S \subseteq [k]} (2\alpha - 1)^{|S|} \hat{w}(S)^2$. The second moment of $X$ is given by
\begin{equation}
\label{eq:E[Xsq]}
\text{E}[X^2] = 2^n \sum_{j = 0}^n {n \choose j} g_w(j/n)^{rn}
\end{equation}
\end{lemma}

The function $g_w(\alpha)$ is similar to the noise sensitivity of a boolean function \cite{o2003computational} and measures the correlation in the value of $w$ between two assignments $\sigma$, $\tau$ which overlap at $\alpha(\sigma, \tau) n$ locations. As a visual example, Figure \ref{fig:g_w} shows how $g_w(\alpha)$ changes for $k$-XORSAT with varying $k$. The key of our proof is showing that $\text{E}[w(\sigma_{I, s})w(\tau_{I, s})] = g_w(\alpha(\sigma, \tau))$ for a random constraint $f_{I, s}$.  

We will now write $\text{E}[X]^2$ in terms of $g_w$. Since $g_w(1/2) = \hat{w}(\emptyset)^2$, plugging this into \eqref{eq:E[X]sq} gives us
\begin{equation}
\label{eq:e[X]sqing}
\text{E}[X]^2 = 2^{2n} g_w(1/2)^{rn} 
\end{equation}
This motivates us to apply the following lemma from \cite{achlioptas2004threshold}, which will allow us to translate bounds on $g_w(\alpha)$ into bounds on $\text{E}[X]^2/\text{E}[X^2]$: 
\begin{lemma}
\label{lem:achlioptaslem}
Let $\phi$ be any real, positive, twice-differentiable function on $[0, 1]$ and let 
\begin{align*}
S_n = \sum_{j = 0}^{n} {n \choose j} \phi(j/n)^n
\end{align*}
Define $\psi$ on $[0, 1]$ as $\psi(\alpha) = \frac{\phi(\alpha)}{\alpha^{\alpha}(1 - \alpha)^{1 - \alpha}}$.
If there exists $\alpha_{\max} \in (0, 1)$ such that $\psi(\alpha_{\max}) > \psi(\alpha)$ for all $\alpha \ne \alpha_{\max}$, and $\psi''(\alpha_{\max}) < 0$, then there exist constants $B, C > 0$ such that for all sufficiently large $n$, 
\begin{equation}
\label{eq:lemachlioptasres}
B \psi(\alpha_{\max})^n \le S_n \le C \psi(\alpha_{\max})^n 
\end{equation}
\end{lemma}
To apply the lemma, we can define $\phi_r(\alpha) = g_w(\alpha)^r$ and $\psi_r(\alpha) = \frac{\phi_r(\alpha)}{\alpha^{\alpha}(1 - \alpha)^{1 - \alpha}}$. Then from \eqref{eq:e[X]sqing}, we note that 
\begin{equation}
\psi_r(1/2)^n = 2^n(g_w(1/2)^r)^n = \text{E}[X]^2/2^n
\label{eq:psiE[X]sq}
\end{equation}
On the other hand, from \eqref{eq:E[Xsq]},
\begin{equation}
\sum_{j = 0}^n {n \choose j} \phi_r(j/n)^n = \text{E}[X^2]/2^n
\label{eq:psiE[Xsq]}
\end{equation}
so if the conditions of Lemma \ref{lem:achlioptaslem} hold for $\alpha_{\max} = 1/2$, we recover that $\text{E}[X]^2/\text{E}[X^2] \ge C$ for some constant $C > 0$. 

One requirement for $\psi_r(\alpha)$ to be maximized at $\alpha = 1/2$ is that $\psi_r'(1/2) = 0$. Expanding $\psi_r'(1/2)$ gives $2g_w(1/2)^{r - 1}(rg_w'(1/2)) = 0$.
Since $g_w(1/2) = \hat{w}(\emptyset)^2 > 0$, we thus require 
\begin{equation}
\label{eq:derivhalf0}
g_w'(1/2) = 2\sum_{S \subseteq [k] : |S| = 1} \hat{w}(S)^2 = 0
\end{equation}
In order to satisfy \eqref{eq:derivhalf0}, we need $\hat{w}(S) = 0$ for all $S \subseteq [k]$ where $|S| = 1$. To use Lemma \ref{lem:achlioptaslem}, we would like to choose $w$ such that \eqref{eq:derivhalf0} holds. We discuss how to choose $w$ to optimize our lower bounds in the appendix. In the next section, we will provide $r$ so that the conditions of Lemma \ref{lem:achlioptaslem} hold at $\alpha = 1/2$ for arbitrary $w$ when \eqref{eq:derivhalf0} is satisfied. 

\subsection{Bounding the Second Moment For Fixed $w$}
We give a general bound on $r$ in terms of our weight function $w$ so that the conditions of Lemma \ref{lem:achlioptaslem} are satisfied for $\alpha = 1/2$. For now, the only constraint we place on $w$ is that \eqref{eq:derivhalf0} holds. The next lemma lets us consider only $\alpha \in [1/2, 1]$.  
\begin{lemma}
\label{lem:alphagehalf}
Let $\alpha \ge 1/2$. Then $g_w(\alpha) \ge g_w(1 - \alpha)$. 
\end{lemma}
This lemma follows because $g_w(\alpha)$ is a polynomial in $(2\alpha - 1)$ with nonnegative coefficients, and $(2\alpha - 1) > 0$ for $\alpha > 1/2$.  

Now we can bound $\psi_r(\alpha)$ for $\alpha \in [1/2, 1]$. Combined with Lemma \ref{lem:alphagehalf}, the next lemma will give conditions on $r$ such that $\psi_r(1/2) > \psi_r(\alpha)$ for all $\alpha \in [0, 1]$. 


\begin{lemma}
\label{lem:rbound}
Let the weight function $w$ satisfy \eqref{eq:derivhalf0}. If 
\begin{equation}
\label{eq:rcond}
r \le \frac{1}{2}\frac{\hat{w}(\emptyset)^2}{\sum_{S : |S| \ge 2} \hat{w}(S)^2}
\end{equation}
$\psi_r(1/2) > \psi_r(\alpha)$ for $\alpha \in [0, 1]$ and $\psi_r''(1/2) < 0$.   
\end{lemma}
Figure \ref{fig:psi} shows how $r$ controls the shape of the function $\psi_r(\alpha)$. As $r$ increases, $\psi_r''(1/2)$ becomes positive and $\psi_r(\alpha)$ will no longer attain a local maximum in that region. The key step in proving Lemma \ref{lem:rbound} is rearranging $\psi_r(1/2) > \psi_r(\alpha)$ and simplify calculations by using approximations for the logarithmic terms that appear. 

Our bound on $r$ compares the average of $w$ over $\{-1, 1\}^k$ with the correlations between $w$ and the Fourier basis functions. If $w$ has strong correlations with the other Fourier basis functions, two assignments which are equal at $\alpha n$ variables will likely either be both satisfying or both not satisfying as $\alpha$ approaches $1$. This increases $\text{E}[X^2]$ but not $\text{E}[X]^2$ and makes Lemma \ref{lem:secondmomentmethod} provide a trivial bound if $r$ is too large. Thus, if $w$ has strong correlations with the Fourier basis functions, we must choose smaller $r$ as reflected by \eqref{eq:rcond}. 

To get the tightest bounds, we wish to maximize the expression in \eqref{eq:rcond}. Although we prove our lemma for general $w$ requiring only \eqref{eq:derivhalf0}, we also need $w(u) = 0$ whenever $f(u) = 0$ to apply our lemma to satisfiability. Recalling our definition of $X$ in \eqref{eq:X}, this condition ensures that $C_f(n, rn)$ has a solution whenever $X \ne 0$. Thus,  
\begin{equation}
\label{eq:wlambdaf}
w(u) = \lambda(u) f(u)
\end{equation} 
for some $\lambda : \{-1, 1\}^k \rightarrow \mathbb{R}$. If we disregard \eqref{eq:derivhalf0}, choosing $\lambda(u) = 1$ would maximize the bound on $r$ in \eqref{eq:rcond}. The additional requirement of \eqref{eq:derivhalf0} for the second moment method to succeed can be viewed as a ``symmetrization penalty'' on $r$. In the appendix, we discuss how to choose $w$ to optimize our bound on $r$ while satisfying \eqref{eq:derivhalf0} and \eqref{eq:wlambdaf}.

\begin{figure*}
\centering
\subfloat[$a=2,b=4$]{
\includegraphics[scale=0.32]{tribes_num_r_20_n_30_num_groups_4_num_trials_50_group_size_2}}
\subfloat[$a=2,b=5$]{
\includegraphics[scale=0.32]{tribes_num_r_20_n_20_num_groups_5_num_trials_50_group_size_2}}
\caption{Proportion of CSPs satisfiable out of 50 trials vs. $r$ for tribes functions. We show our bounds for reference.}
\label{fig:experiments}
\end{figure*}

\subsection{Proving the Main Theorem}
\label{sec:mainproof}
We will combine our lemmas to prove Theorem \ref{thm:main}. 
\begin{proof} [Proof of Theorem \ref{thm:main}] 
We wish to apply the second moment method on $X$ defined in \eqref{eq:X}, where $w$ is a function we use to weigh assignments to individual constraints. We choose $w$ as described in the full version of the paper, which satisfies both \eqref{eq:derivhalf0} and \eqref{eq:wlambdaf}. Since $w$ satisfies \eqref{eq:wlambdaf}, $\Pr[C_f(n, rn) \ \text{is satisfiable}] \ge \Pr[X \ne 0] \ge \text{E}[X]^2/\text{E}[X^2]$
by Lemma \ref{lem:secondmomentmethod}. Now we will use Lemma \ref{lem:rbound} to show that the conditions for Lemma \ref{lem:achlioptaslem} are satisfied for $\phi_r = g_w(\alpha)^r$ and $r$ satisfying \eqref{eq:rcond}.
For our choice of $w$, it follows from the  derivations in the appendix that the RHS of \eqref{eq:rcond} becomes 
\begin{align*}
r < r_{f,\text{low}} = \frac{1}{2} \frac{c}{1 - c} \ \text{where} \ c = \hat{f}(\emptyset) - \frac{\boldsymbol{1}^T A^+ A \boldsymbol{1}}{2^k}
\end{align*}
where $A$ is defined in Section \ref{sec:mainresult}. There is a slight technicality in directly applying Lemma \ref{lem:achlioptaslem} because $\phi_r$ might not be nonnegative for $\alpha < 1/2$; we discuss this in the appendix. Now using \eqref{eq:psiE[X]sq} and \eqref{eq:psiE[Xsq]}, and applying Lemma \ref{lem:achlioptaslem}, we can conclude that there exists $C > 0$ such that 
\begin{align*}
\Pr[C_f(n, rn) \ \text{is satisfiable}] \ge \text{E}[X]^2/\text{E}[X^2] \ge C 
\end{align*}
for sufficiently large $n$ and all $r < r_{f,\text{low}}$.
\end{proof}
There remains a question of what $r_{f,\text{low}}$ we can hope achieve using a second moment method proof where $X$ is defined as in \eqref{eq:X}. The following lemma provides some intuition for this: 
\begin{lemma}
\label{lem:maxr}
In order for the conditions of Lemma \ref{lem:achlioptaslem} to hold at $\alpha_{\max} = 1/2$ for $X$ in the form of \eqref{eq:X} and any choice of $w$ satisfying \eqref{eq:wlambdaf}, we require 
\begin{equation}
\label{eq:maxr}
r < \log 2 / \log \frac{1}{\hat{f}(\emptyset) - \frac{\boldsymbol{1}^T A^+ A \boldsymbol{1}}{2^k}} \le r_{f,\text{up}} = \frac{\log 2}{\log \frac{1}{\hat{f}(\emptyset)}}
\end{equation}
\end{lemma}
 
The difference of $\boldsymbol{1}^T A^+ A \boldsymbol{1}/2^k$ in the lower logarithm compared to $r_{f,\text{up}}$ can be viewed as a ``symmetrization penalty'' necessary for our proof to work. 
While Lemma \ref{lem:maxr} does not preclude applications of the second moment method that do not rely on Lemma \ref{lem:achlioptaslem}, consider what happens for $r$ that do not satisfy \eqref{eq:maxr}. For these $r$, the function $\psi_r(\alpha)$ must obtain a maximum at some $\alpha^* \in [0, 1], \alpha^* \ne 1/2$. If it also happens that $\phi_r(\alpha)$ is nonnegative and twice differentiable on $[0, 1]$, and $\psi_r''(\alpha^*) < 0$, then conditions of Lemma \ref{lem:achlioptaslem} hold, and applying it to $\alpha_{\max} = \alpha^*$ along with \eqref{eq:lemachlioptasres}, \eqref{eq:psiE[X]sq}, and \eqref{eq:psiE[Xsq]} will actually imply that 
\begin{align*}
\frac{\text{E}[X]^2}{\text{E}[X^2]} \le \frac{1}{B} \left(\frac{\psi_r(1/2)}{ \psi_r(\alpha^*)}\right)^n
\end{align*}
for some constant $B > 0$, which gives us an exponentially decreasing, and therefore trivial lower bound for the second moment method. Therefore, we believe that \eqref{eq:maxr} is near the best lower bound on $r_{f,\text{unsat}}$ that we can achieve by applying the second moment method on $X$ in the form of \eqref{eq:X}. 

\section{Experimental Verification of Bounds}
\label{sec:experiments}
We empirically test our bounds with the goal of examining their tightness. For our constraint functions, we will use tribes functions. The tribes function takes the disjunction of $b$ groups of $a$ variables and evaluates to 1 or 0 based on whether the following formula is true: 
\begin{align*}
\text{TRIBES}_{a, b}(x_1, \ldots, x_{ab}) = \vee_{i = 0}^{b - 1} \left(\wedge_{j = 1}^a x_{ia + j}\right)
\end{align*}
where $+1$ denotes true and $-1$ denotes false. For our experiments, we randomly generate CSP formulas based on $\text{TRIBES}_{a, b}$. We use the Dimetheus\footnote{https://www.gableske.net/dimetheus} random CSP solver to solve these formulas, or report if no solution exists. We show our results in Figure \ref{fig:experiments}. As expected, our values for lower bounds $r_{f,\text{low}}$ are looser than the upper bounds $r_{f,\text{up}}$. 

\section{Conclusion}
Using Fourier analysis and the second moment method, we have shown general bounds on $m/n$, the ratio of constraints to variables; for $m/n$ below these bounds, there is constant probability that a random CSP is satisfiable. We demonstrate that our bounds are easily instantiated and can be applied to obtain novel estimates of the satisfiability threshold for many classes of CSPs. Our bounds depend on how easy it is to symmetrize solutions to the constraint function. We provide a heuristic argument to approximate the best possible lower bounds that our application of the second moment method can achieve; these bounds differ from upper bounds on the satisfiability threshold by a ``symmetrization penalty.'' Thus, an interesting direction of future research is to determine whether we can provide tighter upper bounds that account for symmetrization, or whether symmetrization terms are an artificial product of the second moment method.  

\section{Acknowledgements}
This work was supported by the Future of Life Institute
(grant 2016-158687) and by the National Science Foundation (grant 1649208).

\bibliographystyle{aaai}
