Solver for finding fixpoints of a boolean system

boolean-algebrainequalitysystems of equations

Intro: The Problem

My problem relates to solving a system of equations that find the fixpoint of a studied boolean system ($F(X)=X$).

A Simple Example

Let $\bar{x}=\{x_1,x_2,x_3\} \in \{0,1\}$ be some boolean variables of interest.
Let $F=\{f_1,f_2,f_3\}$ be the update functions for these variables (i.e. $x_i(t+1) = f_i(x(t))$), defined as the follow condition functions (they are based on the inputs of each variable – node in the corresponding graph conceptualization):

$f_1 = \begin{cases}
1, -x_2 \ge 0 \\
0, \text{otherwise}
\end{cases}$,
$f_2 = \begin{cases}
1, x_1-x_3 \ge 0 \\
0, \text{otherwise}
\end{cases}$,
$f_3 = \begin{cases}
1, x_1+x_3 \ge 0 \\
0, \text{otherwise}
\end{cases}$

These functions are inspired by this question.

The goal is to to find answer sets $\bar{x}$ (could be zero, one or many per se) for which $F(\bar{x})=\bar{x}$.

The above example is of course a very simple case.
In the end I would like to solve such system of equations with hundreds of variables.
Note that the condition functions will always be linear combinations of each variable's inputs and the variables always boolean.

The Question

I need an efficient solver for this kind of problem (which is known to be NP-hard btw!). E.g. can this problem formulated as constraint programming and solved using Answer Set Programming (ASP) techniques?

Best Answer

So it turns out that ASP can be used to solve this problem! Here I provide a possible encoding of the problem in a file named fp.lp:

% variables                                                                     
var(1..3).                                                                      
                                                                                
% functions                                                                     
% f(Function,Coefficient,Variable)                                              
f(1,-1,2).                                                                      
                                                                                
f(2,1,1).                                                                      
f(2,-1,3).                                                                      
                                                                                
f(3,1,1).                                                                       
f(3,1,3).                                                                       
                                                                                
% guess assignment to variables                                                 
{ init(V) : var(V) }.                                                           
                                                                                
% compute functions                                                             
next(F) :- var(F), #sum { C,V : f(F,C,V), init(V) } >= 0.                       
                                                                                
% check if fixed point                                                          
:- init(V), not next(V).                                                        
:- next(V), not init(V).                                                        
                                                                                
#show next/1.                                                                   
#show init/1.

Using clingo (version 5.4.0) from the command line: clingo fp.lp we get UNSATISFIABLE for this particular instance.

Commenting the fact %f(2,1,1). and running again the clingo solver we get:

Answer: 1
next(3) init(1) init(3) next(1)
SATISFIABLE

Variables that are returned both in the init/1 and next/1 predicates are translated to active boolean variables (1) and those that are missing to inactive values (0). The returned result is thus the boolean vector $\bar{x}=\{x_1=1,x_2=0,x_3=1\}$.

Credits for this answer go to Roland Kaminski. Various other members of the ASP potassco community provided helpful comments and solutions. For more info on ASP, check: https://potassco.org/

Related Solutions

[Math] Functionally complete sets of boolean functions

This is a simple application of Post’s characterization of functional completeness.

Here is an example in three variables below.

Let ${\cal F}=\lbrace f_1,f_2,f_3,f_4 \rbrace$ with

\begin{align} f_1&=0,\\ f_2&=1,\\ f_3(x_1,x_2,x_3)&=x_1\oplus x_2 \oplus x_3,\\ f_4(x_1,x_2,x_3)&=x_1 \wedge x_2 \wedge x_3 \end{align}

Then :

$\lbrace f_1,f_2,f_3 \rbrace$ generates only affine boolean functions, so this set is not functionally complete.

$\lbrace f_1,f_2,f_4 \rbrace$ generates only monotone boolean functions, so this set is not functionally complete.

$\lbrace f_1,f_3,f_4 \rbrace$ generates only falsity-preserving boolean functions, so this set is not functionally complete.

$\lbrace f_2,f_3,f_4 \rbrace$ generates only truth-preserving boolean functions, so this set is not functionally complete.

On the other hand, $\cal F$ is functionally complete (indeed, $\neg x_1$ can be expressed as $f_3(x_1,f_1,f_2)$ and $x_1 \wedge x_2$ as $f_4(f_2,x_1,x_2)$), so it is an irreducible functionally complete set with four functions.

[Math] How to find solutions for four polynomial equations with four unknown variables using RESULTANT THEORY

Short answer: Yes; have a look at, e.g., Algorithms in Real Algebraic Geometry by Basu, Pollack and Roy or Fundamental problems in Algorithmic Algebra by Yap, where the CAD approach for quantifier elimination and U-resultants are introduced.

(Slightly) longer answer, with (hopefully) some intuition on why it works: When computing the resultant w.r.t. $y$ ($x$) of a system of two equations in two variables, you eliminate one variable and conceptually project the solutions in $\mathbb C^2$ onto their $x$- ($y$-)values, which are the zeros of the respective resultant polynomial. You then construct a grid of candidate solutions from those coordinates, and check which points are actually proper solutions.
The same idea allows you to eliminate any variable from any pair of two equations in any number $n$ of variables, and the solutions of the $n-1$-variate elimination polynomial determine the candidate set of solutions for those $n-1$ variables. Do this for all pairs $(f_1, f_2), \dots, (f_{n-1}, f_n)$, and you end up with $n-1$ polynomials in $n-1$ variables. Work you way down until you end up with only one variable, and determine the candidates for this projection; then proceed with another way down to get the next variable.
As in the bivariate case, once you have all the candidates for each coordinate, set up the ($n$-dimensional) grid and verify for each point in the grid whether it is a solution of the original system or not.

To your comment: For simplicity, I sketch the three-variate case over an algebraically closed base field $k$ (e.g., $k=\mathbb{C}$) and assume generic position (that is, no solutions at infinity and that the input polynomials are pairwise coprime).

The resultant $R_{i,j} = \operatorname{res}(f_i,f_j; x_3) \in k[x_1,x_2]$ vanishes at $(\alpha,\beta)$ if and only if there exists an $\gamma \in k$ such that $f_i(\alpha,\beta,\gamma) = f_j(\alpha,\beta,\gamma) = 0$. Hence, the roots of the resultant are the projections of the common zeros of $f_i$ and $f_j$ onto the $x_1x_2$-plane. Those are, in general, different values for $R_{1,2}$, $R_{1,3}$ and $R_{2,3}$.
One step further, $R_{(i,j),(i',j')} = \operatorname{res}(R_{i,j},R_{i',j'}; x_2)$ vanishes at $\alpha$ iff there is a $\beta$ such that $R_{i,j}(\alpha,\beta) = R_{i',j'}(\alpha,\beta) = 0$. This in turn means that there are $\beta, \gamma$ and $\gamma'$ such that $f_i(\alpha,\beta,\gamma) = f_j(\alpha,\beta,\gamma) = 0$ and $f_{i'}(\alpha,\beta,\gamma') = f_{j'}(\alpha,\beta,\gamma') = 0$, but $\gamma$ does not need to be the same value as $\gamma'$! Hence, the $R_{(i,j),(i',j')}$ will only share some roots.

However, the projections of the common solutions of the entire system $f_1 = f_2 = f_3 = 0$ are necessarily a part of all resultant polynomials occuring throughout the procedure. In particular, if you descend down to the one-variate projections onto the $x_i$-axes, the gcd of all the univariate elimination polynomials w.r.t. $x_i$ will vanish at each projection of a solution onto the respective coordinate.

Below you can find some proof-of-concept Sage code which shows the idea with very basic means:

n = 3 # nr of variables (and equations)
d = 3 # degrees of the polynomials
n, d = 2, 4
n, d = 4, [2,2,1,1]

if not type(d) is list:
    d = [d]*n

prec = 1024   # precision used for the verification of solutions (exact computation is too expensive without additional machinery)
CIF = ComplexIntervalField (prec)

R = PolynomialRing (ZZ, 'x', n)
x = R.gens()
f = [ R.random_element (x=-256, y=256, degree=d[i], terms=binomial (n+d[i],n)) for i in range(n) ]

print ("Input:\n================================================================")
for i in range(n):
    print "f{} = {}".format (i, f[i])

import itertools
from sage.rings.polynomial.real_roots import real_roots
from sage.rings.polynomial.complex_roots import complex_roots

isolator = lambda p: real_roots (p, retval='algebraic_real')
isolator = lambda p: p.roots (RealIntervalField (prec)) # much faster, but without guarantees
# isolator = lambda p: complex_roots (p, retval='algebraic_complex') # too slow
isolator = lambda p: p.roots (CIF) # much faster, but without guarantees

res = { (i, ()): f[i] for i in range(n) }

for elim_order in itertools.permutations (range(n)):
    print ("Computing resultant polynomials for elimination order {}...".format (elim_order))
    eliminated = []
    pairs = list (itertools.combinations (range(n), 2))
    for k in elim_order[:-1]:
        eliminated.append (k)
        for pair in pairs:
            p = res[pair[0], tuple(eliminated[:-1])]
            q = res[pair[1], tuple(eliminated[:-1])]
            res[pair, tuple(eliminated)] = p.resultant (q, x[k])
        pairs = list (itertools.combinations (pairs, 2))

elim = {}
for i in range(n):
    print ("Computing gcd of univariate elimination polynomials for variable %s..." % x[i])
    elim[i] = gcd (v for k,v in res.items()
                   if len (k[1]) == n-1 # all but one variable eliminated
                   and i not in k[1])   # x_i left

candidates = {}
for i in range(n):
    print ("Computing roots of the elimination polynomial for %s..." % x[i])
    candidates[i] = isolator (elim[i].univariate_polynomial())

print ("Computing the values of the system's polynomials at the candidate solutions...")
sols = []
for cand in itertools.product (*candidates.itervalues()):
    coord = map (lambda x: x[0], cand)
    coord_CIF = map (CIF, coord)
    mult = map (lambda x: x[1], cand)
    values_CIF = [ f[i](coord_CIF) for i in range(n) ]
    print ("\tF{} ~= {}". format (map (ComplexIntervalField(32), coord_CIF),
                                  map (ComplexIntervalField(32), values_CIF)))
    if all (v.contains_zero() for v in values_CIF):
        print ("\t\t=> report as solution")
        sols.append ((coord, mult))

print ("\n(Approximately verified) solutions:\n================================================================")
for sol in sols:
    print ("{}; mult {}".format (map (ComplexIntervalField(32), sol[0]), sol[1]))