Towards End-to-End ASP Computation

Abstract

We propose an end-to-end approach for Answer Set Programming (ASP) and linear algebraically compute stable models satisfying given constraints. The idea is to implement Lin-Zhao’s theorem together with constraints directly in vector spaces as numerical minimization of a cost function constructed from a matricized normal logic program, loop formulas in Lin-Zhao’s theorem and constraints, thereby no use of symbolic ASP or SAT solvers involved in our approach. We also propose precomputation that shrinks the program size and heuristics for loop formulas to reduce computational difficulty. We empirically test our approach with programming examples including the three-coloring and Hamiltonian cycle problems.

1. Introduction

Computing stable model semantics (Gelfond & Lifshcitz, 1988) lies at the heart of Answer Set Programming (ASP) (Lifschitz, 2008; Marek & Truszczyński, 1999; Niemelä, 1999) and there have been a variety of approaches proposed so far. Early approaches such as smodels (Niemelä & Simons, 1997) used backtracking. Then the concept of loop formula was introduced and approaches that use a Boolean satisfiability problem (SAT) solver to compute stable models based on Lin-Zhao’s theorem (Lin & Zhao, 2004) were proposed. They include ASSAT (Lin & Zhao, 2004) and cmodels (Lierler, 2005) for example. Later more elaborated approaches such as clasp (Gebser et al., 2007, 2012) based on conflict-driven no good learning have been developed. While these symbolic approaches continue to predominate in ASP, there has been another trend towards differentiable methods. For example, differentiable ASP/SAT (Nickles, 2018) computes stable models by an ASP solver that utilizes derivatives of a cost function. More recently, NeurASP (Yang et al., 2020) and SLASH (Skryagin et al., 2022, 2023) combined deep learning and ASP. In their approaches, deep learning is not used in an end-to-end way to compute stable models, but used as a component to compute and learn probabilities represented by special atoms interfacing to ASP. Takemura and Inoue (2024) proposed a neurosymbolic learning pipeline that leverages differentiable computation of supported models. Although their method does not specifically address stable model computation, it bypasses the need for a symbolic solver and illustrates how differentiable computation facilitates integration with deep learning. A step towards end-to-end computation was taken by Aspis et al. (2020) and Takemura and Inoue (2022). They formulated the computation of supported models, a super class of stable models, entirely as fixedpoint computation in vector spaces, and obtain supported models represented by binary vectors. However, there still remains a gap between computing supported models and computing stable models.

In this article, we propose an end-to-end approach for ASP and compute stable models satisfying given constraints in vector spaces. The idea is simple; we implement Lin-Zhao’s theorem (Lin & Zhao, 2004) together with constraints directly in vector spaces as a cost minimization problem, thereby no use of symbolic ASP or SAT solvers involved. Since our approach is numerical and relies solely on vector and matrix operations, future work could explore the potential benefits of parallel computing technologies such as many-core CPUs and GPUs.

Technically, Lin-Zhao’s theorem (Lin & Zhao, 2004) states that a stable model of a ground normal logic program coincides with a supported model which satisfies “loop formulas” associated with the program. Loop formulas are propositional formulas indicating how to get out of infinite loops of top-down rule invocation. We formulate finding such a model as root finding in a vector space of a non-negative cost function represented in terms of the matricized program and loop formulas. The problem is that in whatever approach we may take, symbolic or non-symbolic, computing supported models is NP-hard (for example graph coloring is solved by computing supported models) and there can be exponentially many loop formulas to be satisfied (Lifschitz & Razborov, 2006). We reduce this computational difficulty in two ways. One is precomputation that removes atoms from the search space which are known to be false in any stable model and yields a smaller program. The other is to heuristically choose loop formulas to be satisfied. The latter would mean allowing non-stable model computation, and in our continuous approach, we modify the cost function to be affected only by these chosen loop formulas. The intuition behind this heuristic is that the modified cost function would assign higher cost to models that do not satisfy the chosen loop formulas, thus driving the search process away from them.

Our end-to-end computing framework differs from those by Aspis et al. (2020) and Takemura and Inoue (2022) in that they basically compute supported models and the computing process itself has no mechanism such as loop formulas to exclude non-stable models. In addition, any propositional normal logic program is acceptable in our framework, since we impose no restrictions on the syntax of programs like the multiple definition (MD) condition (Aspis et al., 2020) or the singly defined (SD) condition (Takemura & Inoue, 2022). More importantly, we incorporate the use of constraints, that is, rules with an empty head, which make ASP smooth and practical.

Hence, our contributions include as follows:

–

a proposal of end-to-end computing of stable models in vector spaces for propositional normal logic programs;

–

augmentation of the above by constraints;

–

introduction of precomputation and heuristics to reduce computational difficulty of stable model computation.

We add that since our primary purpose in this article is to establish theoretical feasibility of end-to-end ASP computing in vector spaces, programming examples are small and implementation is of preliminary nature. Furthermore, the main search algorithm we propose in this article is incomplete, in the sense that it does not guarantee reaching a global minimum if it exists, nor it cannot conclusively prove that no solution exists.

In what follows, after preliminaries in Section 2, we formulate the computation of supported models in vector spaces in Section 3 and that of stable models in Section 4. We then show programming examples in Section 5 including ASP programs for the three-coloring problem and the Hamiltonian cycle problem. We there compare performance of precomputation and loop formula heuristics. Section 6 contains related work and Section 7 is the conclusion.

2. Preliminaries

In this article, we mean by a program a propositional normal logic program $P$ which is a finite set of rules of the form $a\leftarrow G$ , where $a$ is an atom called the head, $G$ is a conjunction of literals called the body of the rule, respectively. We equate propositional variables with atoms. A literal is an atom (positive literal) or its negation (negative literal). The logical connective $\mathrm{\neg }$ in this article denotes negation as failure. We suppose $P$ is written in a given set of atoms $\mathcal{A}$ but usually assume $\mathcal{A}$ = atom( $P$ ), that is, the set of atoms occurring in $P$ . We use $G^{+}$ and $G^{-}$ to denote the conjunction of positive and negative literals in $G$ , respectively. $G$ may be empty. The empty conjunction is always true. We call $a\leftarrow G$ rule for $a$ . A rule with an empty head is called a constraint. Let $a\leftarrow G_1,\dots ,a\leftarrow G_m$ be rules for $a$ in $P$ . When $m>0$ , put $\mathrm{iff}(a)=a\iff G_1\vee \cdots \vee G_m$ . When $m=0$ , that is, there is no rule for $a$ , put $\mathrm{iff}(a)=a\iff \mathrm{\perp }$ , where $\mathrm{\perp }$ is a special symbol representing the empty disjunction which is always false. We call $\mathrm{iff}(a)$ the completed rule for $a$ . The completion of $P$ , comp( $P$ ), is defined as comp( $P$ ) = $\{\mathrm{iff}(a)\mid \text{atom }a\text{ occurs in }P\}$ . For a finite set $S$ , we denote the number of elements in $S$ by $|S|$ . So $|P|$ is the number of rules in the program $P$ .

An interpretation (assignment) $I$ over a set of atoms $\mathcal{A}$ is a mapping which determines the truth value of each atom $a\in \mathcal{A}$ . Then the truth value of a formula $F$ is inductively defined by $I$ , and if $F$ becomes true evaluated by $I$ , we say $I$ satisfies $F$ , $F$ is true in $I$ , or $I$ is a model of $F$ and write $I\models F$ . This notation is extended to a set $F=\{F_1,\dots ,F_u\}$ by considering $F$ as a conjunction $F_1\wedge \cdots \wedge F_u$ . For convenience, we always equate $I$ with $\{a\in \mathcal{A}\mid I\models a\}$ , that is, the set of atoms true in $I$ . When $I$ satisfies all rules in the program $P$ , that is, $I\models P$ , $I$ is said to be a model of $P$ . If no rule body contains negative literals, $P$ is said to be a definite program. In that case, $P$ always has the least model (in the sense of set inclusion) $\{a\in \mathcal{A}\mid P⊢a\}$ , that is, the set of atoms provable from $P$ .

A model $I$ of comp( $P$ ) is a supported model of $P$ (Apt et al., 1988; Marek & Subrahmanian, 1992). When $P$ is a definite program, there is at least one supported model, and its least model is also a supported model. In general, there can be multiple supported models for both definite and non-definite programs $P$ . Stable models are a subclass of supported models. They are defined as follows. Given a program $P$ and a model $I$ , remove all rules from $P$ whose body contains a negative literal false in $I$ , then remove all negative literals from the remaining rules. The resulting program, $P^I$ , is called the Gelfond-Lifschitz (GL) reduct of $P$ by $I$ or just the reduct of $P$ by $I$ . It is a definite program and has the least model. If this least model is identical to $I$ , $I$ is called a stable model of $P$ (Gelfond & Lifshcitz, 1988). $P$ may have zero or multiple stable models, as in the case of supported models. Since the existence of a stable model is NP-complete (Marek & Truszczyński, 1999) and so is a supported model, their computation is expected to be hard. Supported models and stable models of a propositional normal logic program coincide when the program is tight (no infinite call chain through positive goals) (Erdem & Lifschitz, 2003; Fages, 1994).

Let $F=d_1\vee \cdots \vee d_h$ be a Boolean formula in $n$ variables (atoms) in disjunctive normal form (DNF), where each $d_i$ ( $1\le i\le h$ ) is a conjunction of literals and called a disjunct of $F$ . When $F$ has no disjunct, $F$ is false.

A walk in a directed graph is a sequence $v_1\to v_2\to \cdots \to v_u$ ( $u\ge 1$ ) of vertices representing the corresponding non-zero sequence of edges $(v_1,v_2),\dots ,(v_{u-1},v_u)$ . When $v_u=v_1$ , it is said to be closed. A cycle is a closed walk $v_1\to v_2\to \cdots \to v_u\to v_1$ , where $\{v_1,\dots ,v_u\}$ are all distinct. A Hamiltonian cycle (HC) is a cycle which visits every vertex exactly once. A path is a walk with no vertex repeated. A directed subgraph is called strongly connected if there are paths from $v_1$ to $v_2$ and from $v_2$ to $v_1$ for any pair of distinct vertices $v_1$ and $v_2$ . This “strongly connected” relation induces an equivalence relation over the set of vertices and an induced equivalence class is called a strongly connected component (SCC).

The positive dependency graph $\mathrm{pdg}(P)$ for a program $P$ is a directed graph where vertices are atoms occurring in $P$ and there is an edge $(a,b)$ from atom $a$ to atom $b$ if and only if (iff) there is a rule $a\leftarrow G$ in $P$ such that $b$ is a positive literal in $G$ . $P$ is said to be tight (Erdem & Lifschitz, 2003; Fages, 1994) ¹ when $\mathrm{pdg}(P)$ is acyclic, that is, has no cycle. A loop $S=\{a_1,\dots ,a_u\}$ $(u>0)$ in $P$ is a set of atoms where for any pair of atoms $a_1$ and $a_2$ in $S$ ( $a_1=a_2$ allowed), there is a path in $\mathrm{pdg}(P)$ from $a_1$ to $a_2$ and also from $a_2$ to $a_1$ . A singleton loop $S=\{a\}$ is induced by a self-referencing rule of the form $a\leftarrow a\wedge G$ where $G$ is possibly empty, that is, a self-loop $a\leftarrow a$ . A support rule for $a$ relative to a loop $S$ is a rule $a\leftarrow G$ such that $G^{+}\cap S=\mathrm{\varnothing }$ . Given a loop $L=\{a_1,\dots ,a_u\}$ and its external support rules $\{a_1\leftarrow G_{11},\dots ,a_1\leftarrow G_{1n},\dots ,a_u\leftarrow G_{u1},\dots ,a_u\leftarrow G_{un}\}$ , the (conjunctive) loop formula is the following implication: $(a_1\wedge \cdots \wedge a_u)\to (G_{11}\vee \cdots \vee G_{1n}\vee \cdots \vee G_{u1}\vee \cdots \vee G_{un})$ .

We denote vectors by bold lower case letters such as $\mathbf{a}$ , where $\mathbf{a}(i)$ represents the $i$ -th element of $\mathbf{a}$ . Vectors are column vectors by default. We use $(\mathbf{a}\bullet \mathbf{b})$ to stand for the inner product (dot product) of vectors $\mathbf{a}$ and $\mathbf{b}$ of the same dimension. $\Vert \mathbf{a}{\Vert }_1$ and $\Vert \mathbf{a}{\Vert }_2$ , respectively, denote the 1-norm and 2-norm of $\mathbf{a}$ , where $\Vert \mathbf{a}{\Vert }_1=\sum |\mathbf{a}(i)|$ and $\Vert \mathbf{a}{\Vert }_2=\sqrt{\sum \mathbf{a}{(i)}^2}$ . We use $\mathbf{1}$ to denote an all-ones vector of appropriate dimension. An interpretation $I$ over a set $\mathcal{A}=\{a_1,\dots ,a_n\}$ of $n$ ordered atoms is equated with an $n$ -dimensional binary vector ${\mathbf{s}}_I\in \{0,1{\}}^n$ such that ${\mathbf{s}}_I(i)=1$ if $a_i$ is true in $I$ and ${\mathbf{s}}_I(i)=0$ otherwise ( $1\le i\le n$ ). ${\mathbf{s}}_I$ is called the vectorized $I$ .

Bold upper case letters such as $\mathbf{A}$ stand for a matrix. We use $\mathbf{A}(i,j)$ to denote the $i,j$ -th element of $\mathbf{A}$ , $\mathbf{A}(i,:)$ the $i$ -th row of $\mathbf{A}$ and $\mathbf{A}(:,j)$ the $j$ -th column of $\mathbf{A}$ , respectively. We often consider one dimensional matrices as (row or column) vectors. $\Vert \mathbf{A}{\Vert }_F$ denotes the Frobenius norm of $\mathbf{A}$ , that is, $\sqrt{\sum _{i,j}\mathbf{A}(i,j{)}^2}$ . Let $\mathbf{A},\mathbf{B}\in {\mathbb{R}}^{m\times n}$ be $m\times n$ matrices. $\mathbf{A}\odot \mathbf{B}$ denotes their Hadamard product, that is, $(\mathbf{A}\odot \mathbf{B})(i,j)=\mathbf{A}(i,j)\mathbf{B}(i,j)$ for $i,j(1\le i\le m,1\le j\le n)$ . $[\mathbf{A};\mathbf{B}]$ designates the $2m\times n$ matrix of $\mathbf{A}$ stacked onto $B$ . We implicitly assume that all dimensions of vectors and matrices in various expressions are compatible. We introduce a piece-wise linear function ${min}_1(x)=min(x,1)$ that returns the lesser of 1 and $x$ as an activation function which is related to the popular activation function $\text{ReLU}(x)=max(x,0)$ by $1-{min}_1(x)=\text{ReLU}(1-x)$ . ${min}_1(\mathbf{A})$ denotes the result of component-wise application of ${min}_1(x)$ to matrix $\mathbf{A}$ . We also introduce thresholding notation. Suppose $\theta$ is a real number and $\mathbf{a}$ an $n$ -dimensional vector. Then $[\mathbf{a}\le \theta ]$ denotes a binary vector obtained by thresholding $\mathbf{a}$ at $\theta$ where for $i(1\le i\le n)$ , $[\mathbf{a}\le \theta ](i)=1$ if $\mathbf{a}(i)\le \theta$ and $[\mathbf{a}\le \theta ](i)=0$ otherwise. $[\mathbf{a}\ge \theta ]$ is treated similarly. We extend thresholding to matrices. Thus $[\mathbf{A}\le 1]$ means a matrix such that $[\mathbf{A}\le 1](i,j)=1$ if $\mathbf{A}(i,j)\le 1$ and $[\mathbf{A}\le 1](i,j)=0$ otherwise. For convenience, we generalize bit inversion to an $n$ -dimensional vector $\mathbf{a}$ and use an expression $\mathbf{1}-\mathbf{a}$ to denote the $n$ -dimensional vector such that $(\mathbf{1}-\mathbf{a})(i)=1-\mathbf{a}(i)$ for $i(1\le i\le n)$ . $\mathbf{1}-\mathbf{A}$ is treated similarly.

3. Computing Supported Models in Vector Spaces

In this section, we formulate the semantics of supported models in vector spaces and show how to compute it by cost minimization.

3.1. Matricized Programs

Definition 1 Matricized program

A program $P$ that has $m$ rules in $n$ atoms is numerically encoded as a pair $\mathbf{P}=(\mathbf{C},\mathbf{D})$ of binary matrices $\mathbf{C}\in \{0,1{\}}^{m\times 2n}$ and $\mathbf{D}\in \{0,1{\}}^{n\times m}$ , which we call a matricized program $P$ .

$\mathbf{C}$ represents rule bodies in $P$ . Suppose atoms are ordered like $\mathcal{A}=\{a_1,\dots ,a_n\}$ and similarly rules are ordered like $\{r_1:a_{i_1}\leftarrow G_1,\dots ,r_m:a_{i_m}\leftarrow G_m\}$ . Then the $j$ -th row $\mathbf{C}(j,:)$ ( $1\le j\le m$ ) encodes the $j$ -th conjunction $G_j$ of the $j$ -th rule $a_{i_j}\leftarrow G_j$ . Write $G_j=a_{i_1}\wedge \cdots \wedge a_{i_p}\wedge \mathrm{\neg }{a}_{i_{p+1}}\wedge \cdots \wedge \mathrm{\neg }{a}_{i_{p+q}}$ ( $1\le p,q\le n$ ). Then an element of $\mathbf{C}(j,:)$ is zero except for $\mathbf{C}(j,i_1)=\cdots =\mathbf{C}(j,i_p)=\mathbf{C}(j,n+i_{p+1})=\cdots =\mathbf{C}(j,n+i_{p+q})=1$ . $\mathbf{D}$ combines these conjunctions as a disjunction (DNF) for each atom in $\mathcal{A}$ . If the $i$ -th atom $a_i\in \mathcal{A}$ ( $1\le i\le n$ ) has rules $\{a_i\leftarrow G_{j_1},\dots ,a_i\leftarrow G_{j_s}\}$ in $P$ , we put $\mathbf{D}(i,j_1)=\cdots =\mathbf{D}(i,j_s)=1$ to represent a disjunction $G_{j_1}\vee \cdots \vee G_{j_s}$ which is the right hand side of the completed rule for $a_i$ : $\mathrm{iff}(a_i)=a_i\iff G_{j_1}\vee \cdots \vee G_{j_s}$ . If $a_i$ has no rule, we put $\mathbf{D}(i,j)=0$ for all $j$ ( $1\le j\le m$ ). Thus the matricized $\mathbf{P}=(\mathbf{C},\mathbf{D})$ can represent the completed program $\text{comp}(P)$ .

For concreteness, we explain by an example below.

Example 1 Encoding a program

Suppose we are given a program $P_0$ below containing three rules $\{r_1,r_2,r_3\}$ in a set of atoms ${\mathcal{A}}_0=\{p,q,r\}$ .

$$\begin{array}{rl}{P}_0& =\begin{array}{l}\{\begin{array}{lllll}p& \leftarrow & q\wedge \mathrm{\neg }r& & \text{: rule }{r}_1\text{ for }p\\ p& \leftarrow & \mathrm{\neg }q& & \text{: rule }{r}_2\text{ for }p\\ q& & & & \text{: rule }{r}_3\text{ for }q\end{array}\end{array}\end{array}$$

(1)

Assuming atoms are ordered as $p,q,r$ and correspondingly so are the rules $\{r_1,r_2,r_3\}$ as in (1), we encode $P_0$ as a pair of matrices ${\mathbf{P}}_0=({\mathbf{C}}_0,{\mathbf{D}}_0)$ . Here ${\mathbf{C}}_0$ represents conjunctions (the bodies of $\{r_1,r_2,r_3\}$ ) and ${\mathbf{D}}_0$ their disjunctions so that they jointly represent $P_0$ .

$$\begin{array}{rl}{\mathbf{C}}_0& =\begin{array}{l}\begin{array}{cccccc}p& q& r& \mathrm{\neg }p& \mathrm{\neg }q& \mathrm{\neg }r\end{array}\\ \left[\begin{array}{cccccc}0& 1& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& \begin{array}{c}\text{: }{r}_1\text{ has the body }q\wedge \mathrm{\neg }r\\ \text{: }{r}_2\text{ has the body }\mathrm{\neg }q\\ \text{: }{r}_3\text{ has the empty body}\end{array}\end{array}\end{array}$$

(2)

$$\begin{array}{rl}{\mathbf{D}}_0& =\begin{array}{l}\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\\ \left[\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]& \begin{array}{c}\text{: }p\text{ has two rules }{r}_1\text{ and }{r}_2\\ \text{: }q\text{ has one rule }{r}_3\\ \text{: }r\text{ has no rule}\end{array}\end{array}\end{array}$$

(3)

As can be seen, ${\mathbf{C}}_0$ represents conjunctions in $P_0$ in such a way that ${\mathbf{C}}_0(1,:)$ , for example, represents the conjunction $q\wedge \mathrm{\neg }r$ of the first rule in $P_0$ by setting ${\mathbf{C}}_0(1,2)={\mathbf{C}}_0(1,6)=1$ , and so on. ${\mathbf{D}}_0$ represents the disjunctions of rule bodies. So ${\mathbf{D}}_0(1,1)={\mathbf{D}}_0(1,2)=1$ means the first atom $p$ in $\{p,q,r\}$ has two rules, the first rule $r_1$ and the second rule $r_2$ , representing a disjunction $(q\wedge \mathrm{\neg }r)\vee \mathrm{\neg }q$ for $p$ .

3.2. Evaluation of Formulas and the Reduct of a Program in Vector Spaces

Here we explain how the propositional formulas and the reduct of a program are evaluated by a model in vector spaces. Let $I$ be a model over a set $\mathcal{A}$ of atoms. Recall that $I$ is equated with a subset of $\mathcal{A}$ . We inductively define the relation “a formula $F$ is true in $I$ ,” $I\models F$ in notation, as follows. For an atom $a$ , $I\models a$ iff $a\in I$ . For a compound formula $F$ , $I\models \mathrm{\neg }F$ iff $I⊭F$ . When $F$ is a disjunction $F_1\vee \cdots \vee F_u$ ( $u\ge 0$ ), $I\models F$ iff there is some $i$ ( $1\le i\le u$ ) s.t. $I\models F_i$ . So the empty disjunction ( $u=0$ ) is always false. We consider a conjunction $F_1\wedge \cdots \wedge F_u$ as a syntax sugar for $\mathrm{\neg }(\mathrm{\neg }{F}_1\vee \cdots \vee \mathrm{\neg }{F}_u)$ using De Morgan’s law. Consequently, the empty conjunction is always true. Let $P$ be a program having $m$ ordered rules in $n$ ordered atoms as before and $G=a_{i_1}\wedge \cdots \wedge a_{i_p}\wedge \mathrm{\neg }{a}_{i_{p+1}}\wedge \cdots \wedge \mathrm{\neg }{a}_{i_{p+q}}$ the body of a rule $a\leftarrow G$ in $P$ . By definition, $I\models G$ ( $G$ is true in $I$ ) iff $\{a_{i_1},\dots ,a_{i_p}\}\subseteq I$ and $\{a_{i_{p+1}},\dots ,a_{i_{p+q}}\}\cap I=\mathrm{\varnothing }$ . Also let $\mathrm{iff}(a_i)=a_i\iff G_{j_1}\vee \cdots \vee G_{j_s}$ be the completed rule for an atom $a_i$ in $P$ . We see $I\models \mathrm{iff}(a_i)$ iff $(a_i\in I\text{iff}I\models G_{j_1}\vee \cdots \vee G_{j_s})$ .

Now we isomorphically embed the above symbolic evaluation to the one in vector spaces. Let $I$ be a model over ordered atoms $\mathcal{A}=\{a_1,\dots ,a_n\}$ . We first vectorize $I$ as a binary column vector ${\mathbf{s}}_I$ such that ${\mathbf{s}}_I(i)=1$ if $a_i\in I$ and ${\mathbf{s}}_I(i)=0$ ( $1\le i\le n$ ) otherwise, and introduce the dualized ${\mathbf{s}}_I$ written as ${\mathbf{s}}_I^{\delta }$ by ${\mathbf{s}}_I^{\delta }=[{\mathbf{s}}_I;(\mathbf{1}-{\mathbf{s}}_I)]$ . ${\mathbf{s}}_I^{\delta }$ is a vertical concatenation of ${\mathbf{s}}_I$ and the bit inversion of ${\mathbf{s}}_I$ .

Consider a matricized program $\mathbf{P}=(\mathbf{C},\mathbf{D})$ $(\mathbf{C}\in \{0,1{\}}^{m\times 2n}$ , $\mathbf{D}\in \{0,1{\}}^{n\times m})$ and its $j$ -th rule $r_j$ having a body $G_j$ represented by $\mathbf{C}(j,:)$ . Compute $\mathbf{C}(j,:){\mathbf{s}}_I^{\delta }$ which is the number of true literals in $I$ in $G_j$ and compare it with the number of literals $|\mathbf{C}(j,:){|}_1$ ² in $G_j$ . When $|\mathbf{C}(j,:){|}_1=\mathbf{C}(j,:){\mathbf{s}}_I^{\delta }$ holds, all literals in $G_j$ are true in $I$ and hence the body $G_j$ is true in $I$ . In this way, we can algebraically compute the truth value of each rule body, but since we consider a conjunction as a negated disjunction, we instead compute $\mathbf{C}(j,:)(\mathbf{1}-{\mathbf{s}}_I^{\delta })$ which is the number of false literals in $G_j$ . If this number is non-zero, $G_j$ have at least one literal false in $I$ , and hence $G_j$ is false in $I$ . The converse is also true. The existence of a false literal in $G_j$ is thus computed by ${min}_1(\mathbf{C}(j,:)(\mathbf{1}-{\mathbf{s}}_I^{\delta }))$ which is $1$ if there is a false literal, and $0$ otherwise. Consequently, $1-{min}_1(\mathbf{C}(j,:)(\mathbf{1}-{\mathbf{s}}_I^{\delta }))=1$ if there is no false literal in $G_j$ and vice versa. In other words, $1-{min}_1(\mathbf{C}(j,:)(\mathbf{1}-{\mathbf{s}}_I^{\delta }))$ computes $I\models G_j$ .

Now let $\{a_i\leftarrow G_{j_1},\dots ,a_i\leftarrow G_{j_s}\}$ be an enumeration of rules for $a_i\in \mathcal{A}$ and $G_{j_1}\vee \cdots \vee G_{j_s}$ the disjunction of the rule bodies. $d_i=\sum _{t=1}^s(1-{min}_1(\mathbf{C}(j_t,:)(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ is the number of rule bodies in $\{G_{j_1},\dots ,G_{j_s}\}$ that are true in $I$ . Noting $\mathbf{D}(i,j)=1$ if $j\in \{j_1,\dots ,j_s\}$ and $\mathbf{D}(i,j)=0$ otherwise by construction of $\mathbf{D}$ in $\mathbf{P}=(\mathbf{C},\mathbf{D})$ , we replace the summation $\sum _{t=1}^s$ by matrix multiplication and obtain $d_i=\mathbf{D}(i,:)(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ . Introduce a column vector ${\mathbf{d}}_I=\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ . We have ${\mathbf{d}}_I(i)=d_i$ = the number of rules for $a_i$ whose bodies are true in $I$ ( $1\le i\le n$ ).

In the case of ${\mathbf{P}}_0=({\mathbf{C}}_0,{\mathbf{D}}_0)$ in (1) having three rules $\{r_1,r_2,r_3\}$ , take a model $I_0=\{p,q\}$ over the ordered atom set ${\mathcal{A}}_0=\{p,q,r\}$ , where $p$ and $q$ are true in $I_0$ but $r$ is false in $I_0$ . Then we have ${\mathbf{s}}_{I_0}=[110{]}^T$ , ${\mathbf{s}}_{I_0}^{\delta }=[110001{]}^T$ , $\mathbf{1}-{\mathbf{s}}_{I_0}^{\delta }=[001110{]}^T$ , and finally, ${\mathbf{C}}_0(\mathbf{1}-{\mathbf{s}}_{I_0}^{\delta })=[010{]}^T$ . The last equation says that the rule bodies of $r_1$ , $r_2$ , and $r_3$ have, respectively, zero, one, and zero literal false in $I_0$ . Hence ${min}_1({\mathbf{C}}_0(\mathbf{1}-{\mathbf{s}}_{I_0}^{\delta }))=[010{]}^T$ indicates that only the second rule body is false and the other two bodies are true in $I_0$ . So its bit inversion $\mathbf{1}-{min}_1({\mathbf{C}}_0(\mathbf{1}-{\mathbf{s}}_{I_0}^{\delta }))=[101{]}^T$ indicates that the second rule body is false in $I_0$ , while others are true in $I_0$ . Thus by combining these truth values in terms of disjunctions $D_0$ , we obtain ${\mathbf{d}}_{I_0}={\mathbf{D}}_0(\mathbf{1}-{min}_1({\mathbf{C}}_0(\mathbf{1}-{\mathbf{s}}_{I_0}^{\delta })))=[110{]}^T$ . The elements in ${\mathbf{d}}_{I_0}=[110{]}^T$ denote for each atom $a\in {\mathcal{A}}_0$ the number of rules for $a$ whose body is true in $I_0$ . For example, ${\mathbf{d}}_{I_0}(1)=1$ means that the first atom $p$ in ${\mathcal{A}}_0$ has one rule ( $p\leftarrow q\wedge \mathrm{\neg }r$ ) whose body ( $q\wedge \mathrm{\neg }r$ ) is true in $I_0$ . Likewise, ${\mathbf{d}}_{I_0}(2)=1$ means that the second atom $q$ has one rule ( $q\leftarrow$ ) whose body (empty) is true in $I_0$ . ${\mathbf{d}}_{I_0}(3)=0$ indicates that the third atom $r$ has no such rule. Therefore, ${min}_1({\mathbf{d}}_{I_0})=[110{]}^T$ denotes the truth values of the right hand sides of the completed rules $\{\mathrm{iff}(p),\mathrm{iff}(q),\mathrm{iff}(r)\}$ evaluated by $I_0$ .

Proposition 1

Let $\mathbf{P}=(\mathbf{C},\mathbf{D})$ be a matricized program $P$ in a set of atoms $\mathcal{A}$ and ${\mathbf{s}}_I$ a vectorized model $I$ over $\mathcal{A}$ . Put ${\mathbf{d}}_I=\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ . It holds that

$$I\models \text{comp}(P)\text{iff}\Vert {\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I){\Vert }_2=0.$$

(4)

Proof.

Put $n=|\mathcal{A}|$ . Suppose $I\models \text{comp}(P)$ and write $\mathrm{iff}(a_i)$ , the completed rule for an atom $a_i\in \mathcal{A}$ ( $1\le i\le n$ ), as $\mathrm{iff}(a_i)=a_i\iff G_{j_1}\vee \cdots \vee G_{j_s}$ ( $s\ge 0$ ). We have $I\models \mathrm{iff}(a_i)$ . So if ${\mathbf{s}}_I(i)=1$ , $I\models a_i$ , and hence $I\models G_{j_1}\vee \cdots \vee G_{j_s}$ , giving $d_i\ge 1$ because $d_i$ is the number of rule bodies in $\{G_{j_1},\dots ,G_{j_s}\}$ that are true in $I$ . So ${min}_1(d_i)=1$ holds. Otherwise, if ${\mathbf{s}}_I(i)=0$ , we have $I⊭a_i$ and $I⊭G_{j_1}\vee \cdots \vee G_{j_s}$ . Consequently, none of the rule bodies are true in $I$ and we have $d_i={min}_1(d_i)=0$ . Putting the two together, we have ${\mathbf{s}}_I(i)=d_i$ . Since $i$ is arbitrary, we conclude ${\mathbf{s}}_I={\min }_1({\mathbf{d}}_I)$ , or equivalently $\Vert {\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I){\Vert }_2=0$ . The converse is similarly proved.

Proposition 1 says that whether $I$ is a supported model of the program $P$ or not is determined by computing ${\mathbf{s}}_I-{min}_1({\mathbf{d}}_I)$ in vector spaces whose complexity is $O(mn)$ where $m$ is the number of rules in $P$ , $n$ that of atoms occurring in $P$ . In the case of ${\mathbf{P}}_0=({\mathbf{C}}_0,{\mathbf{D}}_0)$ with ${\mathbf{s}}_{I_0}=[110{]}^T$ as in the aforementioned example, since ${\mathbf{s}}_{I_0}={min}_1({\mathbf{d}}_{I_0})$ holds, it follows from Proposition 1 that $I_0$ is a supported model of $P_0$ .

We next show how $P^I$ , the reduct of $P$ by $I$ , is dealt with in vector spaces. We assume $P$ has $m$ rules $\{r_1,\dots ,r_m\}$ with a set $\mathcal{A}=\{a_1,\dots ,a_n\}$ of $n$ ordered atoms as before. We first show the evaluation of the reduct of the matricized program $\mathbf{P}=(\mathbf{C},\mathbf{D})$ by a vectorized model ${\mathbf{s}}_I$ . Write $\mathbf{C}\in \{0,1{\}}^{m\times 2n}$ as $\mathbf{C}=[{\mathbf{C}}^{pos}{\mathbf{C}}^{neg}]$ , where ${\mathbf{C}}^{pos}\in \{0,1{\}}^{m\times n}$ (resp. ${\mathbf{C}}^{neg}\in \{0,1{\}}^{m\times n}$ ) is the left half (resp. the right half) of $\mathbf{C}$ representing the positive literals (resp. negative literals) of each rule body in $\mathbf{C}$ . Compute ${\mathbf{M}}^{neg}=\mathbf{1}-{min}_1({\mathbf{C}}^{neg}{\mathbf{s}}_I)$ . It is an $m\times 1$ matrix (treated as a column vector here) such that ${\mathbf{M}}^{neg}(j)=0$ if the body of $r_j$ contains a negative literal false in $I$ and ${\mathbf{M}}^{neg}(j)=1$ otherwise ( $1\le j\le m$ ). Let $r_j^{+}$ be a rule $r_j$ with negative literals in the body deleted. We see that $P^I=\{r_j^{+}\mid {\mathbf{M}}^{neg}(j)=1,1\le j\le m\}$ and $P^I$ is syntactically represented by $({\mathbf{C}}^{pos},{\mathbf{D}}^I)$ , where ${\mathbf{D}}^I=\mathbf{D}$ with columns $\mathbf{D}(:,j)$ replaced by the zero column vector if ${\mathbf{M}}^{neg}(j)=0$ ( $1\le j\le m$ ). ${\mathbf{D}}^I(i,:)$ denotes a rule set $\{r_j^{+}\mid {\mathbf{D}}^I(i,j)=1,1\le j\le m\}$ in $P^I$ for $a_i\in \mathcal{A}$ . We call $P^I=({\mathbf{C}}^{pos},{\mathbf{D}}^I)$ the matricized reduct of $P$ by $I$ .

The matricized reduct $P^I=({\mathbf{C}}^{pos},{\mathbf{D}}^I)$ is evaluated in vector spaces as follows. Compute ${\mathbf{M}}^{pos}={\mathbf{M}}^{neg}\odot (\mathbf{1}-{min}_1({\mathbf{C}}^{pos}(\mathbf{1}-{\mathbf{s}}_I)))$ . ${\mathbf{M}}^{pos}$ denotes the truth values of rule bodies in $P^I$ evaluated by $I$ . Thus ${\mathbf{M}}^{pos}(j)=1$ ( $1\le j\le m$ ) if $r_j^{+}$ is contained in $P^I$ and its body is true in $I$ . Otherwise ${\mathbf{M}}^{pos}(j)=0$ and $r_j^{+}$ is not contained in $P^I$ or the body of $r_j^{+}$ is false in $I$ . Introduce ${\mathbf{d}}_I^{+}=\mathbf{D}{\mathbf{M}}^{pos}$ . ${\mathbf{d}}_I^{+}(i)$ ( $1\le i\le n$ ) is the number of rules in $P^I$ for the $i$ -th atom $a_i$ in $\mathcal{A}$ whose bodies are true in $I$ .

Proposition 2

Let $\mathbf{P}=(\mathbf{C},\mathbf{D})$ be a matricized program $P$ in a set $\mathcal{A}=\{a_1,\dots ,a_n\}$ of $n$ ordered atoms and $I$ a model over $\mathcal{A}$ . Write $\mathbf{C}=[{\mathbf{C}}^{pos}{\mathbf{C}}^{neg}]$ as above. Let ${\mathbf{s}}_I$ be the vectorized model $I$ . Compute ${\mathbf{M}}^{neg}=\mathbf{1}-{\min }_1({\mathbf{C}}^{neg}{\mathbf{s}}_I)$ , ${\mathbf{M}}^{pos}={\mathbf{M}}^{neg}\odot (\mathbf{1}-{\min }_1({\mathbf{C}}^{pos}(\mathbf{1}-{\mathbf{s}}_I)))$ , and ${\mathbf{d}}_I^{+}=D{M}^{pos}$ . Also compute ${\mathbf{d}}_I=\mathbf{D}(\mathbf{1}-{\min }_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ . Then, $I\models \mathrm{comp}(P)$ , $\Vert {\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I){\Vert }_2=0$ , $\Vert {\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I^{+}){\Vert }_2=0$ , and $I\models \mathrm{comp}(P^I)$ are all equivalent (proof in Appendix A.1).

From the viewpoint of end-to-end ASP, Proposition 2 means that we can obtain a supported model $I$ as a binary solution ${\mathbf{s}}_I$ of the equation ${\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I)=0$ derived from $P$ or ${\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I^{+})=0$ derived from the reduct $P^I$ . Either equation is possible and gives the same result but their computation will be different. This is because the former equation ${\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I)$ is piecewise linear w.r.t. ${\mathbf{s}}_I$ , whereas the latter ${\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I^{+})$ is piecewise quadratic w.r.t. ${\mathbf{s}}_I$ .

Example 2 Evaluation of a reduct

Now look at $P_0=\{r_1,r_2,r_3\}$ in (1) and a model $I_0=\{p,q\}$ again. $\begin{array}{lcl}{P}_0^{I_0}& =& \begin{array}{l}\{\begin{array}{lll}p& \leftarrow & q\\ q& \leftarrow & \end{array}\end{array}\end{array}$ is the reduct of $P_0$ by $I_0$ . $P_0^{I_0}$ has the least model $\{p,q\}$ that coincides with $I_0$ . So $I_0$ is a stable model of $P_0$ . To simulate the reduction process of $P_0$ in vector spaces, let ${\mathbf{P}}_0=({\mathbf{C}}_0,{\mathbf{D}}_0)$ be the matricized $P_0$ . We first decompose ${\mathbf{C}}_0$ in (2) as ${\mathbf{C}}_0=[{\mathbf{C}}_0^{pos}{\mathbf{C}}_0^{neg}]$ , where ${\mathbf{C}}_0^{pos}$ is the positive part and ${\mathbf{C}}_0^{neg}$ the negative part of ${\mathbf{C}}_0$ . They are

$$\begin{array}{r}{\mathbf{C}}_0^{pos}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\text{and}{\mathbf{C}}_0^{neg}=\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 0& 0& 0\end{array}\right].\end{array}$$

Let ${\mathbf{s}}_{I_0}=[110{]}^T$ be the vectorized $I_0$ . We first compute ${\mathbf{M}}_0^{neg}=\mathbf{1}-{min}_1({\mathbf{C}}_0^{neg}{\mathbf{s}}_{I_0})$ to determine rules to be removed. Since ${\mathbf{M}}_0^{neg}=[101{]}^T$ , the second rule $r_2$ , indicated by ${\mathbf{M}}_0^{neg}(2)=0$ , is removed from $P_0$ , giving $P_0^{I_0}=\{r_1^{+},r_3^{+}\}$ . Using ${\mathbf{M}}_0^{neg}$ and ${\mathbf{D}}_0$ shown in (3), we then compute ${\mathbf{M}}_0^{pos}={\mathbf{M}}_0^{neg}\odot (\mathbf{1}-{min}_1({\mathbf{C}}_0^{pos}(\mathbf{1}-{\mathbf{s}}_{I_0})))=[101{]}^T$ and ${\mathbf{d}}_{I_0}^{+}={\mathbf{D}}_0{\mathbf{M}}_0^{pos}=[110{]}^T$ . The elements of ${\mathbf{d}}_{I_0}^{+}$ denote the number of rule bodies in $P_0^{I_0}$ that are true in $I_0$ for each atom. Thus, since ${\mathbf{s}}_{I_0}={min}_1({\mathbf{d}}_{I_0}^{+})(=[110{]}^T)$ holds, $I_0$ is a supported model of $P_0$ by Proposition 2.

3.3. Cost Minimization for Supported Models

Having linear algebraically reformulated several concepts in logic programming, we tackle the problem computing supported models in vector spaces. Although there already exist approaches for this problem, we tackle it without assuming any condition on programs while allowing constraints. Aspis et al. formulated the problem as solving a non-linear equation containing a sigmoid function (Aspis et al., 2020). They encode normal logic programs differently from ours based on Sakama’s encoding (Sakama et al., 2017) and impose the MD condition on programs which is rather restrictive. No support is provided for constraints in their approach. Later Takemura and Inoue proposed another approach (Takemura & Inoue, 2022) which encodes a program in terms of a single matrix and evaluates conjunctions by the number of true literals. They compute supported models by minimizing a non-negative function, not solving an equation like (Aspis et al., 2020). Their programs are, however, restricted to those satisfying the SD condition and constraints are not considered.

Here we introduce an end-to-end way of computing supported models in vector spaces through cost minimization of a new cost function based on the evaluation of disjunction. We impose no syntactic restriction on programs and allow constraints. We believe that these two features make our end-to-end ASP approach more feasible.

We can base our supported model computation either on Proposition 1 or on Proposition 2. In the latter case, we have to compute GL reduction which requires complicated computation compared to the former case. So for the sake of simplicity, we explain the former. Then our task in vector spaces is to find a binary vector ${\mathbf{s}}_I$ representing a supported model $I$ of a matricized program $\mathbf{P}=(\mathbf{C},\mathbf{D})$ that satisfies $\Vert {\mathbf{s}}_I-{\min }_1({\mathbf{d}}_I){\Vert }_2=0$ , where ${\mathbf{d}}_I=\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}_I^{\delta })))$ . For this task, we relax ${\mathbf{s}}_I\in \{0,1{\}}^n$ to $\mathbf{s}\in {\mathbb{R}}^n$ and introduce a non-negative cost function $L_{SU}$ :

$$L_{SU}=0.5\cdot (\Vert \mathbf{s}-{\min }_1(\mathbf{d}){\Vert }_2^2+{\ell }_2\cdot \Vert \mathbf{s}\odot (\mathbf{1}-\mathbf{s}){\Vert }_2^2),\text{where}{\ell }_2>0\text{and}\mathbf{d}=\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }))).$$

(5)

Proposition 3

Let $L_{SU}$ be defined from a program $\mathbf{P}=(\mathbf{C},\mathbf{D})$ as above.

$L_{SU}=0$ iff $\mathbf{s}$ is a binary vector representing a supported model of $P$ .

Proof.

Apparently if $L_{SU}=0$ , we have $\Vert \mathbf{s}-{\min }_1(\mathbf{d}){\Vert }_2^2=0$ and $\Vert \mathbf{s}\odot (\mathbf{1}-\mathbf{s}){\Vert }_2^2=0$ . The second equation means $\mathbf{s}$ is binary ( $x(1-x)=0\iff x\in \{0,1\}$ ), and the first equation means this binary $\mathbf{s}$ is a vector representing a supported model of $P$ by Proposition 1. The converse is obvious.

$L_{SU}$ is piecewise differentiable and we can obtain a supported model of $P$ as a root $\mathbf{s}$ of $L_{SU}$ by minimizing $L_{SU}$ to zero using Newton’s method. The Jacobian $J_{a_{SU}}$ required for Newton’s method is derived as follows. We assume $P$ is written in $n$ ordered atoms $\{a_1,\dots ,a_n\}$ and $\mathbf{s}=[u_1,\dots ,u_n{]}^T$ represents their continuous truth values, where $\mathbf{s}(p)=s_p\in \mathbb{R}$ is the continuous truth value for atom $a_p$ ( $1\le p\le n$ ). For the convenience of derivation, we introduce the dot product $(\mathbf{A}\bullet \mathbf{B})=\sum _{i,j}\mathbf{A}(i,j)\mathbf{B}(i,j)$ of matrices $\mathbf{A}$ and $\mathbf{B}$ and a one-hot vector ${\mathbf{I}}_p$ which is a zero vector except for the $p$ -th element and ${\mathbf{I}}_p(p)=1$ . We note $(\mathbf{A}\bullet (\mathbf{B}\odot \mathbf{C}))=((\mathbf{B}\odot \mathbf{A})\bullet \mathbf{C})$ and $(\mathbf{A}\bullet (\mathbf{B}\mathbf{C}))=(({\mathbf{B}}^T\mathbf{A})\bullet \mathbf{C})=((\mathbf{A}{\mathbf{C}}^T)\bullet \mathbf{B})$ hold (see Appendix B.1 for details).

Let $\mathbf{P}=(\mathbf{C},\mathbf{D})$ be the matricized program and write $\mathbf{C}=[{\mathbf{C}}^{pos}{\mathbf{C}}^{neg}]$ . Introduce $\mathbf{N}$ , $\mathbf{M}$ , $\mathbf{d}$ , $\mathbf{E}$ , $\mathbf{F}$ , and compute $L_{SU}$ by

$$\begin{array}{r}\begin{array}{lclll}\mathbf{N}& =& \mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta })={\mathbf{C}}^{pos}(\mathbf{1}-\mathbf{s})+{\mathbf{C}}^{neg}\mathbf{s}& & \text{: (continuous) counts of false literals in the rule bodies}\\ \mathbf{M}& =& \mathbf{1}-{\min }_1(\mathbf{N})& & \text{: (continuous) truth values of the rule bodies}\\ \mathbf{d}& =& \mathbf{D}\mathbf{M}& & \text{: (continuous) counts of true disjuncts for each atom}\\ \mathbf{E}& =& {\min }_1(\mathbf{d})-\mathbf{s}& & \text{: error between the estimated truth values of atoms and }\mathbf{s}\\ \mathbf{F}& =& \mathbf{s}\odot (\mathbf{1}-\mathbf{s})& & \text{: (continuous) 0 iif }\mathbf{s}\text{ is binary }\\ L_{sq}& =& (\mathbf{E}\bullet \mathbf{E})\\ L_{nrm}& =& (\mathbf{F}\bullet \mathbf{F})\\ L_{SU}& =& 0.5\cdot (L_{sq}+{\ell }_2\cdot L_{nrm}).\end{array}\end{array}$$

(6)

We then compute the Jacobian $J_{a_{SU}}$ of $L_{SU}$ as follows (full derivation in Appendix B.2):

$$\begin{array}{rl}{J}_{a_{SU}}& =\left(\frac{\mathrm{\partial }{L}_{sq}}{\mathrm{\partial }\mathbf{s}}\right)+{\ell }_2\cdot \left(\frac{\mathrm{\partial }{L}_{nrm}}{\mathrm{\partial }\mathbf{s}}\right)\\ & =({\mathbf{C}}^{pos}-{\mathbf{C}}^{neg}{)}^T([\mathbf{N}\le 1]\odot ({\mathbf{D}}^T([\mathbf{d}\le 1]\odot \mathbf{E})))-\mathbf{E}+{\ell }_2\cdot ((\mathbf{1}-2\mathbf{s})\odot \mathbf{F}),\end{array}$$

(7)

$$\begin{array}{rl}& \text{where}\mathbf{N}=\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }),\mathbf{d}=\mathbf{D}(\mathbf{1}-{\min }_1(\mathbf{N})),\mathbf{E}={\min }_1(\mathbf{d})-\mathbf{s},\text{and}\mathbf{F}=\mathbf{s}\odot (\mathbf{1}-\mathbf{s}).\end{array}$$

Note that since $\mathbf{s}$ is a vector, the Jacobian in this case is also a vector.

3.4. Adding Constraints

A rule which has no head like $\leftarrow a\wedge \mathrm{\neg }b$ is called a constraint. We oftentimes need supported models which satisfy constraints. Since constraints are just rules without a head, we encode constraints as rule bodies in a program using a binary matrix $\hat{\mathbf{C}}=[{\hat{\mathbf{C}}}^{pos}{\hat{\mathbf{C}}}^{neg}]$ . We call $\hat{\mathbf{C}}$ constraint matrix. We introduce ${\mathbf{N}}_{\hat{\mathbf{c}}}$ , a non-negative function $L_{\hat{\mathbf{c}}}$ of $\mathbf{s}$ and $L_{\hat{\mathbf{c}}}$ ’s Jacobian $J_{a_{\hat{\mathbf{c}}}}$ as follows (derivation in Appendix B.3):

$$\begin{array}{rl}{\mathbf{N}}_{\hat{\mathbf{c}}}& =\hat{\mathbf{C}}(\mathbf{1}-{\mathbf{s}}^{\delta })={\hat{\mathbf{C}}}^{pos}(\mathbf{1}-\mathbf{s})+{\hat{\mathbf{C}}}^{neg}\mathbf{s}\text{: number of literals falsified by }\mathbf{s}\\ L_{\hat{\mathbf{c}}}& =(\mathbf{1}\bullet (\mathbf{1}-{\min }_1({\mathbf{N}}_{\hat{\mathbf{c}}}))),\text{where }\mathbf{1}\text{ is an all-ones vector}\text{: counts of violated constraints}\end{array}$$

(8)

$$\begin{array}{rl}{J}_{a_{\hat{\mathbf{c}}}}& =({\hat{\mathbf{C}}}^{pos}-{\hat{\mathbf{C}}}^{neg}{)}^T[{\mathbf{N}}_{\hat{\mathbf{c}}}\le 1]\text{: the Jacobian of }{L}_{\hat{\mathbf{c}}}\end{array}$$

(9)

The meaning of ${\mathbf{N}}_{\hat{\mathbf{c}}}$ and $L_{\hat{\mathbf{c}}}$ is clear when $\mathbf{s}$ is binary. Note that any binary $\mathbf{s}$ is considered as a model over a set $\mathcal{A}=\{a_1,\dots ,a_n\}$ of $n$ ordered atoms in an obvious way. Suppose $k$ constraints are given to be satisfied. Then $\hat{\mathbf{C}}$ is a $k\times 2n$ binary matrix and ${\mathbf{N}}_{\hat{\mathbf{c}}}$ is a $k\times 1$ matrix. ${\mathbf{N}}_{\hat{\mathbf{c}}}(i)$ ( $1\le i\le k$ ) is the number of literals falsified by $\mathbf{s}$ in a conjunction $G_i$ of the $i$ -th constraint $\leftarrow G_i$ . So ${\mathbf{N}}_{\hat{\mathbf{c}}}(i)=0$ , or equivalently $1-{\min }_1({\mathbf{N}}_{\hat{\mathbf{c}}}(i))=1$ implies $G_i$ has no false literal, that is, $\mathbf{s}\models G_i$ , and vice versa. Hence $L_{\hat{\mathbf{c}}}=\sum _{i=1}^k(1-{\min }_1({\mathbf{N}}_{\hat{\mathbf{c}}}(i)))=(\mathbf{1}\bullet (\mathbf{1}-{\min }_1({\mathbf{N}}_{\hat{\mathbf{c}}})))$ equals the number of violated constraints. Consequently, when $\mathbf{s}$ is binary, we can say that $L_{\hat{\mathbf{c}}}=0$ iff all constraints are satisfied by $\mathbf{s}$ .

When $\mathbf{s}$ is not binary but just a real vector $\mathbf{s}\in {\mathbb{R}}^n$ , ${\mathbf{N}}_{\hat{\mathbf{c}}}$ and $L_{\hat{\mathbf{c}}}$ are thought to be a continuous approximation to their binary counterparts. Since $L_{\hat{\mathbf{c}}}$ is a piecewise differentiable non-negative function of $\mathbf{s}$ , the approximation error can be minimized to zero by Newton’s method using $J_{a_{\hat{\mathbf{c}}}}$ in (9).

3.5. An Algorithm for Computing Supported Models With Constraints

Here we present a minimization algorithm for computing supported models of the matricized program $\mathbf{P}=(\mathbf{C},\mathbf{D})$ which satisfy constraints represented by a constraint matrix $\hat{\mathbf{C}}$ . We first combine $L_{SU}$ and $L_{\hat{\mathbf{c}}}$ into $L_{SU+\hat{\mathbf{c}}}$ using ${\ell }_3>0$ .

$$\begin{array}{rl}{L}_{SU+\hat{\mathbf{c}}}& =L_{SU}+{\ell }_3\cdot L_{\hat{\mathbf{c}}}\end{array}$$

(10)

$$\begin{array}{rl}& =0.5\cdot (\Vert \mathbf{s}-{\min }_1(\mathbf{d}){\Vert }_2^2+{\ell }_2\cdot \Vert \mathbf{s}\odot (\mathbf{1}-\mathbf{s}){\Vert }_2^2)+{\ell }_3\cdot (\mathbf{1}\bullet (\mathbf{1}-{\min }_1(\hat{\mathbf{C}}(1-{\mathbf{s}}^{\delta })))){\ell }_2>0,{\ell }_3>0,\\ & \text{where}\mathbf{d}=\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }))),\\ J_{a_{SU+\hat{\mathbf{c}}}}& =J_{a_{SU}}+{\ell }_3\cdot J_{a_{\hat{\mathbf{c}}}}.\end{array}$$

(11)

The next proposition is immediate from Proposition 3.

Proposition 4

$L_{SU+\hat{\mathbf{c}}}=0$ iff   $\mathbf{s}$ represents a supported model of $P$ satisfying a constraint matrix $\hat{\mathbf{C}}$ .

We compute $L_{SU}$ in $L_{SU+\hat{\mathbf{c}}}$ by (6) and $L_{\hat{\mathbf{c}}}$ by (8), and their Jacobians $J_{a_{SU}}$ and $J_{a_{\hat{\mathbf{c}}}}$ by (7) and by (9), respectively. We minimize the non-negative $L_{SU+\hat{\mathbf{c}}}$ to zero by Newton’s method using Algorithm 1. It finds a solution ${\mathbf{s}}_{*}$ of $L_{SU+\hat{\mathbf{c}}}=0$ which represents a supported model of $P$ satisfying constraint matrix $\hat{\mathbf{C}}$ . The updating formula is derived from the first-order Taylor expansion of $L_{SU+\hat{\mathbf{c}}}$ and by solving $L_{SU+\hat{\mathbf{c}}}+(J_{a_{SU*c}}\bullet ({\mathbf{s}}_{\mathrm{new}}-\mathbf{s}))=0$ w.r.t. ${\mathbf{s}}_{\mathrm{new}}$ . The updating formula with a learning rate $\alpha >0$ is thus defined as follows:

$${\mathbf{s}}_{\mathrm{new}}=\mathbf{s}-\alpha \left(\frac{L_{SU+\hat{\mathbf{c}}}}{(J_{a_{SU+\hat{\mathbf{c}}}}\bullet J_{a_{SU+\hat{\mathbf{c}}}})}\right)J_{a_{SU+\hat{\mathbf{c}}}}.$$

(12)

Algorithm 1 is a double-loop algorithm where the inner $j$ -loop updates $\mathbf{s}\in {\mathbb{R}}^n$ repeatedly to minimize $L_{SU+\hat{\mathbf{c}}}$ while thresholding $\mathbf{s}$ into a binary solution candidate ${\mathbf{s}}_{*}\in \{0,1{\}}^n$ for $L_{SU+\hat{\mathbf{c}}}=0$ . The outer $i$ -loop is for retry when the inner loop fails to find a solution. The initialization at line 3 is carried out by sampling $\mathbf{s}(i)\sim \mathcal{N}(0,1)+0.5$ ( $1\le i\le n$ ) where $\mathcal{N}(0,1)$ is the standard normal distribution. Lines 6 to 8 collectively perform thresholding of $\mathbf{s}$ into a binary ${\mathbf{s}}_{*}$ . As the inner loop repeats, $L_{SU+\hat{\mathbf{c}}}$ becomes smaller and smaller and so do $L_{sq}$ and $L_{nrm}$ in $L_{SU}$ . $L_{sq}$ being small means $\mathbf{s}$ is close to a supported model of $P$ while $L_{nrm}$ being small means each element of $\mathbf{s}$ is close to $\{0,1\}$ . So binarization ${\mathbf{s}}_{*}=[\mathbf{s}\ge \theta ]$ with an appropriate threshold $\theta$ ³ has a good chance of yielding a binary ${\mathbf{s}}_{*}$ representing a supported model of $P$ satisfying constraints represented by $\hat{\mathbf{C}}$ . It may happen that the inner loop fails to find a solution. In such a case, we retry another $j$ -loop with perturbated $\mathbf{s}$ at line 12. There $\mathbf{s}$ is perturbated by $\mathbf{s}\leftarrow 0.5(\mathbf{s}+\mathrm{\Delta }+0.5)$ where $\mathrm{\Delta }\sim \mathcal{N}(0,1)$ before the next $j$ -loop.

4. Computing Stable Models in Vector Spaces

4.1. Loop Formulas and Stable Models

Let $\mathbf{P}=(\mathbf{C},\mathbf{D})$ be a matricized program in a set of atoms $\mathcal{A}=\{a_1,\dots ,a_n\}$ having $m$ rules $\{a_{i_1}\leftarrow G_1,\dots ,a_{i_m}\leftarrow G_m\}$ , where $\mathbf{C}\in \{0,1{\}}^{m\times 2n}$ and $\mathbf{D}\in \{0,1{\}}^{n\times m}$ . We assume atoms and rules are ordered as indicated.

Computing a supported model of $P$ is equivalent to computing any binary fixedpoint $\mathbf{s}\in \{0,1{\}}^n$ such that $\mathbf{s}={\min }_1(\mathbf{D}(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }))))$ in vector spaces and in this sense, it is conceptually simple (though NP-hard). Contrastingly since stable models are a proper subclass of supported models, if one wishes to obtain precisely stable models through fixedpoint computation, the exclusion of non-stable models is necessary. Lin-Zhao’s theorem (Lin & Zhao, 2004) states that $I$ is a stable model of $P$ iff $I$ is a supported model of $P$ and satisfies a set of formulas called loop formulas associated with $P$ .

Let $S=\{h_1,\dots ,h_p\}\subseteq \mathcal{A}$ be a loop in $P$ . Recall that $S$ is a set of atoms which are strongly connected in the positive dependency graph of $P$ . ⁴ A support rule for $h$ with respect to $S$ is a rule $h\leftarrow H$ such that $H^{+}\cap S=\mathrm{\varnothing }$ . $H$ is called a support body for $S$ . Introduce a (conjunctive) loop formula for $S$ by

$$LF(S)=(h_1\wedge \cdots \wedge h_p)\to (H_1\vee \cdots \vee H_q),\text{where }\{H_1,\dots ,H_q\}\text{ are support bodies for }S.$$

(13)

Then define loop formulas associated with $P$ as $LF(P)=\{LF(S)\mid S\text{ is a loop in }P\}$ , which is treated as the conjunction of its elements. We note that in the original form (Lin & Zhao, 2004), the antecedent of $LF(S)$ is a disjunction $(h_1\vee \cdots \vee h_p)$ . Later it was shown that the disjunctive and conjunctive loop formulas are equivalent (Ferraris et al., 2006), and we choose to use the conjunctive form of $LF(S)$ as it is easier to satisfy using our method. We evaluate $LF(P)$ by a real vector $\mathbf{s}\in {\mathbb{R}}^n$ . Introduce an external support matrix ${\mathbf{E}}_{sup}\in \{0,1{\}}^{n\times m}$ by ${\mathbf{E}}_{sup}(i,j)=1$ if there is a support rule $a_i\leftarrow G_j$ for $a_i\in \mathcal{A}$ , else ${\mathbf{E}}_{sup}(i,j)=0$ ( $1\le i\le n,1\le j\le m$ ). Suppose there are $w$ loops $\{S_1,\dots ,S_w\}$ in $P$ . Introduce a loop matrix ${\mathbf{L}}_{oop}\in \{0,1{\}}^{w\times m}$ such that ${\mathbf{L}}_{oop}(v,j)=1$ if the $v$ -th loop $S_v$ has $G_j$ as a support body for $S_v$ , else ${\mathbf{L}}_{oop}(v,j)=0$ ( $1\le v\le w$ ).

Example 3 Encoding loop formulas

Suppose we are given a program $P_{L0}$ :

$$\begin{array}{rl}{P}_{L0}& =\begin{array}{l}\{\begin{array}{lllll}p& \leftarrow & q\wedge \mathrm{\neg }r& & \text{: rule }{r}_1\text{ for }p\\ p& \leftarrow & \mathrm{\neg }s& & \text{: rule }{r}_2\text{ for }p\\ q& \leftarrow & p& & \text{: rule }{r}_3\text{ for }q\\ r& \leftarrow & r& & \text{: rule }{r}_4\text{ for }r\end{array}\end{array}\end{array}$$

(14)

This program contains two loops: $S_1=\{p,q\}$ and $S_2=\{r\}$ . In this case, only $S_1$ has an external support body $\mathrm{\neg }s$ . Thus, the external support matrix ${\mathbf{E}}_{sup0}$ and the loop matrix ${\mathbf{L}}_{oop0}$ for this program are as follows:

$$\begin{array}{rl}{\mathbf{E}}_{sup0}& =\begin{array}{l}\begin{array}{cccc}{r}_1& r_2& r_3& r_4\end{array}\\ \left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]& \begin{array}{c}\text{: }p\text{ has the body of }{r}_2\text{ as the support body}\\ \text{: }q\text{ has no external support}\\ \text{: }r\text{ has no external support}\\ \text{: }s\text{ is not part of any loops}\end{array}\end{array}\end{array}$$

(15)

$$\begin{array}{rl}{\mathbf{L}}_{oop0}& =\begin{array}{l}\begin{array}{cccc}{r}_1& r_2& r_3& r_4\end{array}\\ \left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right]& \begin{array}{c}\text{: }{S}_1\text{ has the body of }{r}_2\text{ as its support body}\\ \text{: }{S}_2\text{ has no support bodies}\end{array}\end{array}\end{array}$$

(16)

We then introduce a loss function $L_{LF}$ , which is a non-negative piecewise linear function of $\mathbf{s}$ .

$$\begin{array}{rl}\mathbf{M}& =\mathbf{1}-{\min }_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }))\text{: (continuous) truth values by }\mathbf{s}\text{ of the rule bodies in }P\\ {\mathbf{S}}_v& ={\mathbf{L}}_{oop}(v,:)\text{: represents the }v\text{-th loop in }\{{\mathbf{S}}_1,\dots ,{\mathbf{S}}_w\}\\ A_v& ={\mathbf{S}}_v(\mathbf{1}-\mathbf{s})+{\mathbf{S}}_v{\mathbf{E}}_{sup}\mathbf{M}\text{: (continuous) counts of true disjuncts by }\mathbf{s}\text{ of }LF(S_v)\\ L_{LF}& =\sum _{v=1}^w(1-{\min }_1(A_v))\end{array}$$

(17)

Proposition 5

Let $L_{LF}$ be defined as above. When $\mathbf{s}$ is a binary vector representing a model $I$ over $\mathcal{A}$ , it holds that $L_{LF}=0$ iff $\mathbf{s}\models LF(P)$ .

Proof.

Suppose $L_{LF}=0$ and $\mathbf{s}$ is binary. A summand $(1-{\min }_1(A_v))$ in $L_{LF}$ (17) corresponds to the $v$ -th loop $S_v=\{h_1,\dots ,h_p\}$ and is non-negative. Consider $LF(S_v)=(h_1\wedge \cdots \wedge h_p)\to (H_1\vee \cdots \vee H_q)$ as a disjunction $\mathrm{\neg }{h}_1\vee \cdots \vee \mathrm{\neg }{h}_p\vee H_1\vee \cdots \vee H_q$ . Then $L_{LF}=0$ implies $(1-{\min }_1(A_v))=0$ , or equivalently $A_v\ge 1$ . Consequently, as $\mathbf{s}$ is binary, we have ${\mathbf{S}}_v(\mathbf{1}-\mathbf{s})\ge 1$ or ${\mathbf{S}}_v{\mathbf{E}}_{sup}\mathbf{M}\ge 1$ . The former means $I\models \mathrm{\neg }{h}_1\vee \cdots \vee \mathrm{\neg }{h}_p$ . The latter, ${\mathbf{S}}_v{\mathbf{E}}_{sup}\mathbf{M}\ge 1$ , means $I\models H_1\vee \cdots \vee H_q$ . This is because the element $({\mathbf{E}}_{sup}\mathbf{M})(i)$ is the number of support rules for $a_i\in \mathcal{A}$ whose bodies are true in $I$ $\mathbf{s}$ ( $1\le i\le n$ ), and hence ${\mathbf{S}}_v{\mathbf{E}}_{sup}\mathbf{M}\ge 1$ means some support body $H_r$ ( $1\le r\le q$ ) for $S_v$ is true in $I$ . So in either case $I\models LF(S_v)$ . Since $v$ is arbitrary, we have $I\models LF(P)$ . The converse is straightforward and omitted.

The Jacobian $J_{a_{LF}}$ of $L_{LF}$ is computed as follows (derivation in Appendix B.4):

$$\begin{array}{rl}\mathbf{N}& =\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta }),\\ N_v& ={\mathbf{S}}_v(\mathbf{1}-\mathbf{s}),\\ M_v& ={\min }_1(N_v),\\ J_{a_{LF}}& =\frac{\mathrm{\partial }{L}_{LF}}{\mathrm{\partial }\mathbf{s}}=\sum _{v=1}^w-\left(\frac{\mathrm{\partial }{\min }_1(A_v)}{\mathrm{\partial }\mathbf{s}}\right)\\ & =\sum _{v=1}^w[A_v\le 1]([N_v\le 1]{\mathbf{S}}_v^T+((({\mathbf{S}}_v{\mathbf{E}}_{sup})\odot [\mathbf{N}\le 1{]}^T)({\mathbf{C}}^{neg}-{\mathbf{C}}^{pos}{)}^T).\end{array}$$

(18)

Here $\mathbf{C}=[{\mathbf{C}}^{pos}{\mathbf{C}}^{neg}]$ and ${\mathbf{S}}_v$ , $A_v$ , and $\mathbf{M}$ are computed by (17).

Now introduce a new cost function $L_{SU+\hat{\mathbf{c}}+LF}$ by (19) that incorporates $L_{LF}$ and compute its Jacobian $J_{a_{SU+\hat{\mathbf{c}}+LF}}$ by (11).

$$\begin{array}{rl}{L}_{SU+\hat{\mathbf{c}}+LF}& =L_{SU+\hat{\mathbf{c}}}+{\ell }_4\cdot L_{LF},\text{where }{\ell }_4>0,\end{array}$$

(19)

$$\begin{array}{rl}{J}_{a_{SU+\hat{\mathbf{c}}+LF}}& =J_{a_{SU+\hat{\mathbf{c}}}}+{\ell }_4\cdot J_{a_{LF}}.\end{array}$$

(20)

By combining Propositions 4, 5, and Lin-Zhao’s theorem (Lin & Zhao, 2004), the following is obvious.

Proposition 6

$\mathbf{s}$ is a stable model of $P$ satisfying constraints represented by $\hat{\mathbf{C}}$ iff $\mathbf{s}$ is a root of $L_{SU+\hat{\mathbf{c}}+LF}$ .

We compute such $\mathbf{s}$ by Newton’s method using Algorithm 1 with a modified update rule (12) such that $L_{SU+\hat{\mathbf{c}}}$ and $J_{a_{SU+\hat{\mathbf{c}}}}$ are replaced by $L_{SU+\hat{\mathbf{c}}+LF}$ and $J_{a_{SU+\hat{\mathbf{c}}+LF}}$ , respectively.

When a program $P$ is tight (Fages, 1994), for example, when rules have no positive literal in their bodies, $P$ has no loop and hence $LF$ is empty. In such a case, we directly minimize $L_{SU+\hat{\mathbf{c}}}$ instead of using $L_{SU+\hat{\mathbf{c}}+LF}$ with the empty $LF$ .

4.2. LF Heuristics

Minimizing $L_{SU+\hat{\mathbf{c}}+LF}$ is a general way of computing stable models under constraints. It is applicable to any program and gives us a theoretical framework for computing stable models in an end-to-end way without depending on symbolic systems. However, there can be exponentially many loops and they make the computation of $L_{LF}$ (17) extremely difficult or practically impossible. To mitigate this seemingly insurmountable difficulty, we propose two heuristics which use a subset of loop formulas.

$L{F}_{max}$ :

The first heuristic is $L{F}_{max}$ . We consider only a set $L{F}_{max}$ of loop formulas associated with SCCs in the positive dependency graph $\mathrm{pdg}(P)=(V,E)$ of a program $P$ . In the case of a singleton SCC $\{a\}$ , $a$ must have a self-loop in $\mathrm{pdg}(P)$ . We compute SCCs in $O(|E|+|V|)$ time by Tarjan’s algorithm (Tarjan, 1972).

$L{F}_{min}$ :

In this heuristic, instead of SCCs (maximal strongly connected subgraphs), we choose minimal strongly connected subgraphs, that is, cycle graphs. Denote by $L{F}_{min}$ the set of loop formulas associated with cycle graphs in $\mathrm{pdg}(P)$ . We use an enumeration algorithm described by Liu and Wang (2006) to enumerate cycles and construct $L{F}_{min}$ due to its simplicity.

We remark that although $L{F}_{max}$ and $L{F}_{min}$ can exclude some of non-stable models, they do not necessarily exclude all of non-stable models. However, the role of loop formulas in our framework is entirely different from the one in symbolic ASP. Namely, the role of LF in our framework is not to logically reject non-stable models but to guide the search process by their gradient information in the continuous search space. Hence, we expect, as actually observed in experiments in the next section, some loop formulas have the power of guiding the search process to a root of $L_{SU+\hat{\mathbf{c}}+LF}$ .

4.3. Precomputation

We introduce here precomputation. The idea is to remove atoms from the search space which are false in every stable model. It downsizes the program and realizes faster model computation.

When a program $P$ in a set $\mathcal{A}$ = atom( $P$ ) is given, we transform $P$ to a definite program $P^{+}$ by removing all negative literals from the rule bodies in $P$ . Since $P^{+}\supseteq P^I$ holds as a set of rules for any model $I$ , we have $LM(P^{+})\supseteq LM(P^I)$ , where $LM(P)$ denotes the least model of a definite program $P$ . When $I$ is a stable model, $LM(P^I)=I$ holds and we have $LM(P^{+})\supseteq I$ . By taking the complements of both sides, we can say that if an atom $a$ is outside of $LM(P^{+})$ , that is, if $a$ is false in $LM(P^{+})$ , so is $a$ in any stable model $I$ of $P$ . Thus, by precomputing the least model $LM(P^{+})$ , we can remove a set of atoms ${\mathcal{F}}_P=\mathcal{A}\setminus LM(P^{+})$ from our consideration as they are known to be false in any stable model. We call ${\mathcal{F}}_P$ stable false atoms. Of course, this precomputation needs additional computation of $LM(P^{+})$ but it can be done in linear time proportional to the size of $P^{+}$ , that is, the total number of occurrences of atoms in $P^{+}$ (Dowling & Gallier, 1984). ⁵ Accordingly precomputing the least model $LM(P^{+})$ makes sense if the benefit of removing stable false atoms from the search space outweighs linear time computation for $LM(P^{+})$ , which is likely to happen when we deal with programs with positive literals in the rule bodies.

More concretely, given a program $P$ and a set of constraints $K$ , we can obtain downsized ones, $P^{\prime }$ and $K^{\prime }$ , as follows.

Step 1:

Compute the least model $LM(P^{+})$ and the set of stable false atoms ${\mathcal{F}}_P=\text{atom}(P)\setminus LM(P^{+})$ .

Step 2:

Define

$$\begin{array}{rl}{G}^{\prime }& =\text{conjunction }G\text{ with negative literals }\{\mathrm{\neg }a\in G\mid a\in {\mathcal{F}}_P\}\text{ removed}\\ P^{\prime }& =\{a\leftarrow G^{\prime }\mid a\leftarrow G\in P,a\notin {\mathcal{F}}_P,G^{+}\cap {\mathcal{F}}_P=\mathrm{\varnothing }\},\text{where }{G}^{+}=\text{positive literals in }G,\end{array}$$

(21)

$$\begin{array}{rl}{K}^{\prime }& =\{\leftarrow G^{\prime }\mid \leftarrow G\in K,G^{+}\cap {\mathcal{F}}_P=\mathrm{\varnothing }\}.\end{array}$$

(22)

Proposition 7

Let $P^{\prime }$ and $K^{\prime }$ be, respectively, the program (21) and constraints (22). Also let $I^{\prime }$ be a model over atom( $P^{\prime }$ ). Expand $I^{\prime }$ to a model $I$ over atom( $P$ ) by assuming every atom in ${\mathcal{F}}_P$ is false in $I$ . Then

$I^{\prime }$ is a stable model of $P^{\prime }$ satisfying constraints $K^{\prime }$ iff $I$ is a stable model of $P$ satisfying constraints $K$ .

Proof.

We prove first $I^{\prime }$ is a stable model of $P^{\prime }$ iff $I$ is a stable model of $P$ . To prove it, we prove $LM(P^{\prime I^{\prime }})=LM(P^I)$ as set.

Let $a\leftarrow G^{\prime +}$ be an arbitrary rule in $P^{\prime I^{\prime }}$ . Correspondingly there is a rule $a\leftarrow G^{\prime }$ in $P^{\prime }$ such that $I^{\prime }\models G^{\prime -}$ . So there is a rule $a\leftarrow G$ in $P$ such that $G^{\prime }=G\setminus \{\mathrm{\neg }b\mid b\in {\mathcal{F}}_P\}$ and $G^{+}\cap {\mathcal{F}}_P=\mathrm{\varnothing }$ . $I^{\prime }\models G^{\prime -}$ implies $I\models G^{-}$ by construction of $I$ from $I^{\prime }$ . So $a\leftarrow G^{+}$ is contained in $P^I$ , which means $a\leftarrow G^{\prime +}$ is contained in $P^I$ because $G^{\prime +}=G^{+}$ (recall that $G^{\prime }=G\setminus \{\mathrm{\neg }b\mid b\in {\mathcal{F}}_P\}$ and $G^{\prime }$ and $G$ have the same set of positive literals). Thus since $a\leftarrow G^{\prime +}$ is an arbitrary rule, we conclude $P^{\prime I^{\prime }}\subseteq P^I$ , and hence $LM(P^{\prime I^{\prime }})\subseteq LM(P^I)$ .

Now consider $a\in LM(P^I)$ . There is a selective linear definite clause (SLD) derivation for $\leftarrow a$ from $P^I$ . Let $b\leftarrow G^{+}\in P^I$ be a rule used in the derivation which is derived from the rule $b\leftarrow G\in P$ such that $I\models G^{-}$ . Since $P^I\subseteq P^{+}$ , we have $LM(P^I)\subseteq LM(P^{+})$ and hence $LM(P^I)\cap {\mathcal{F}}_P=\mathrm{\varnothing }$ , that is, $LM(P^I)$ contains no stable false atom. So $b\notin {\mathcal{F}}_P$ and $G^{+}\cap {\mathcal{F}}_P=\mathrm{\varnothing }$ because every atom in the SLD derivation must belong in $LM(P^I)$ . Accordingly $b\leftarrow G^{\prime }\in P^{\prime }$ . On the other hand, $I\models G^{-}$ implies $I^{\prime }\models G^{\prime -}$ . So $b\leftarrow G^{\prime }$ is in $P^{\prime }$ and $b\leftarrow G^{\prime +}$ is in $P^{\prime I^{\prime }}$ . Therefore $b\leftarrow G^{+}$ is in $P^{\prime I^{\prime }}$ because $G^{\prime +}=G^{+}$ . Thus every rule used in the derivation for $\leftarrow a$ from $P^I$ is also a rule contained in $P^{\prime I^{\prime }}$ , which means $a\in LM(P^{\prime I^{\prime }})$ . Since $a$ is arbitrary, it follows that $LM(P^I)\subseteq LM(P^{\prime I^{\prime }})$ . By putting $LM(P^{\prime I^{\prime }})\subseteq LM(P^I)$ and $LM(P^I)\subseteq LM(P^{\prime I^{\prime }})$ together, we conclude $LM(P^I)=LM(P^{\prime I^{\prime }})$ .

Then, if $I^{\prime }$ is a stable model of $P^{\prime }$ , we have $I^{\prime }=LM(P^{\prime I^{\prime }})=LM(P^I)$ as set. Since $I=I^{\prime }$ as set, we have $I=LM(P^I)$ as set, which means $I$ is a stable model of $P$ . Likewise when $I$ is a stable model of $P$ , we have $I=LM(P^I)=LM(P^{\prime I^{\prime }})$ and $I=I^{\prime }$ as set. So $I^{\prime }=LM(P^{\prime I^{\prime }})$ as set and $I^{\prime }$ is a stable model of $P^{\prime }$ .

As for the constraints, consider a constraint $\leftarrow G^{\prime }$ in $K^{\prime }$ . We consider two cases, where $b\in {\mathcal{F}}_P$ occurs positively and negatively in $G$ and $G^{\prime }$ . In the former case, the body remains the same between $G^{\prime }$ and $G$ , thus if $I^{\prime }\models G^{\prime }$ then $I\models G$ and vice versa. In the latter case, because $\mathrm{\neg }b$ always evaluates to true, the negative occurrence $\mathrm{\neg }b$ in the body of the constraint will not change the result of the conjunction $G$ . Thus, combining the two cases and since the constraint $\leftarrow G^{\prime }$ is arbitrary, we conclude that $I^{\prime }\models K^{\prime }$ iff $I\models K$ .

5. Programming Examples

In this section, we apply our ASP approach to examples as a proof of concept and examine the effectiveness of precomputation and heuristics. Since large scale computing is out of scope in this article, the program size is mostly small. ⁶

5.1. The Three-Coloring Problem

We first deal with the three-coloring problem. Suppose we are given a graph $G_1$ . The task is to color the vertices of the graph blue, red, and green so that no two adjacent vertices have the same color like (b) as shown in Figure 1.

Figure 1.

Three-coloring problem: (a) graph $G_1$ and (b) three-coloring.

There are four nodes $\{a,b,c,d\}$ in the graph $G_1$ . We assign a set of three color atoms (Boolean variables) to each node to represent their color. For example, node $a$ is assigned three color atoms $\{a_1(\mathrm{red}),a_2(\mathrm{blue}),a_3(\mathrm{green})\}$ . We need to represent two facts by these atoms.

• –

Each node has a unique color chosen from $\{red,blue,green\}$ . So color atoms assigned to each node are in an XOR relation. We represent this fact by a tight program $P_1$ below containing three rules for each node.

$$\begin{array}{rl}{P}_1& =\{\begin{array}{lll}{a}_1\leftarrow \mathrm{\neg }{a}_2\wedge \mathrm{\neg }{a}_3,& a_2\leftarrow \mathrm{\neg }{a}_3\wedge \mathrm{\neg }{a}_1,& a_3\leftarrow \mathrm{\neg }{a}_1\wedge \mathrm{\neg }{a}_2\\ b_1\leftarrow \mathrm{\neg }{b}_2\wedge \mathrm{\neg }{b}_3,& b_2\leftarrow \mathrm{\neg }{b}_3\wedge \mathrm{\neg }{b}_1,& b_3\leftarrow \mathrm{\neg }{b}_1\wedge \mathrm{\neg }{b}_2\\ c_1\leftarrow \mathrm{\neg }{c}_2\wedge \mathrm{\neg }{c}_3,& c_2\leftarrow \mathrm{\neg }{c}_3\wedge \mathrm{\neg }{c}_1,& c_3\leftarrow \mathrm{\neg }{c}_1\wedge \mathrm{\neg }{c}_2\\ d_1\leftarrow \mathrm{\neg }{d}_2\wedge \mathrm{\neg }{d}_3,& d_2\leftarrow \mathrm{\neg }{d}_3\wedge \mathrm{\neg }{d}_1,& d_3\leftarrow \mathrm{\neg }{d}_1\wedge \mathrm{\neg }{d}_2\end{array}\end{array}$$

(23)

–

Two nodes connected by an edge must have a different color. We represent this fact in terms of constraints.

$$\begin{array}{rl}{K}_1& =\{\begin{array}{lll}\leftarrow a_1\wedge b_1,& \leftarrow a_2\wedge b_2,& \leftarrow a_3\wedge b_3\\ \leftarrow a_1\wedge c_1,& \leftarrow a_2\wedge c_2,& \leftarrow a_3\wedge c_3\\ \leftarrow b_1\wedge c_1,& \leftarrow b_2\wedge c_2,& \leftarrow b_3\wedge c_3\\ \leftarrow b_1\wedge d_1,& \leftarrow b_2\wedge d_2,& \leftarrow b_3\wedge d_3\\ \leftarrow d_1\wedge c_1,& \leftarrow d_2\wedge c_2,& \leftarrow d_3\wedge c_3\end{array}\end{array}$$

(24)

Assuming an ordering of atoms $\{a_1,a_2,a_3,\dots ,d_1,d_2,d_3\}$ , the normal logic program $P_1$ shown in (23) is matricized to ${\mathbf{P}}_1=({\mathbf{C}}_1,{\mathbf{D}}_1)$ , where ${\mathbf{D}}_1$ is a $(12\times 12)$ binary identity matrix (because there are 12 atoms and each atom has just one rule) and ${\mathbf{C}}_1$ is a $(12\times 24)$ binary matrix shown in (25). Constraints listed in (24) are a matricized to a ( $15\times 12$ ) constraint matrix ${\hat{\mathbf{C}}}_{K_1}$ (26). In (25) and (26), $a$ , for example, stands for a triple $(a_1{a}_2{a}_3)$ and $\mathrm{\neg }a$ for $(\mathrm{\neg }{a}_1\mathrm{\neg }{a}_2\mathrm{\neg }{a}_3)$ .

$$\begin{array}{rl}{\mathbf{C}}_1& =\begin{array}{l}\begin{array}{cccccccc}a& b& c& d& \mathrm{\neg }a& \mathrm{\neg }b& \mathrm{\neg }c& \mathrm{\neg }d\end{array}\\ \left[\begin{array}{ccccccccc}& & & & {\mathbf{H}}_3& & & & \\ & & & & & {\mathbf{H}}_3& & & \\ & & & & & & {\mathbf{H}}_3& & \\ & & & & & & & {\mathbf{H}}_3& \end{array}\right],& \text{where}{\mathbf{H}}_3=\left[\begin{array}{ccc}0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right].\end{array}\end{array}$$

(25)

$$\begin{array}{rl}{\hat{\mathbf{C}}}_{K_1}& =\begin{array}{l}\begin{array}{cccc}a& b& c& d\end{array}\\ \left[\begin{array}{cccc}{\mathbf{E}}_3& {\mathbf{E}}_3& & \\ {\mathbf{E}}_3& & {\mathbf{E}}_3& \\ & {\mathbf{E}}_3& {\mathbf{E}}_3& \\ & {\mathbf{E}}_3& & {\mathbf{E}}_3\\ & & {\mathbf{E}}_3& {\mathbf{E}}_3\end{array}\right],& \text{where}{\mathbf{E}}_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].\end{array}\end{array}$$

(26)

We run Algorithm 1 on program $P_1$ with constraints $K_1$ to find a supported model (solution) of $P_1$ satisfying $K_1$ . ⁷

To measure time to find a model, we conduct 10 trials ⁸ of running Algorithm 1 with $max\mathrm{\_}try=20$ , $max\mathrm{\_}itr=50$ , ${\ell }_2={\ell }_3=0.1$ and take the average. The result is 0.104 s (0.070) ⁹ on average. Also to check the ability of finding different solutions, we perform 10 trials of Algorithm 1 ¹⁰ and count the number of different solutions in the returned solutions. #solutions in Table 1 is the average of 10 such measurements. Considering there are six solutions and the naive implementation, the number of different solutions found by the algorithm, which was 5.2 on average, seems rather high.

Table 1.

Time and the Number of Solutions.

Time(s)	# solutions
6.7 (0.7)	5.2 (0.9)

Next we check the scalability of our approach by a simple problem. We consider the three-coloring of a cycle graph like (a) in Figure 2. In general, given a cycle graph that has $n$ nodes, we encode its three-coloring problem as in the previous example by a matricized program ${\mathbf{P}}_2=({\mathbf{C}}_2,{\mathbf{D}}_2)$ and a constraint matrix ${\hat{\mathbf{C}}}_{K_2}$ , where ${\mathbf{D}}_2(3n\times 3n)$ is an identity matrix and $\mathbf{C}(3n\times 6n)$ and ${\hat{\mathbf{C}}}_{K_2}(3n\times 6n)$ represent, respectively, rules and constraints. There are $2^n+2(-1{)}^n$ solutions ( $n\ge 3$ ) in $2^{3n}$ possible assignments for $3n$ atoms. ¹¹ So the problem will be exponentially difficult as $n$ goes up.

Figure 2.

Convergence and scalability: (a) cycle graph, (b) minimization of $L_{SU+\hat{\mathbf{c}}}$ with retry, and (c) scalability.

The graph (b) in Figure 2 is an example of convergence curve of $L_{SU+\hat{\mathbf{c}}}$ by Algorithm 1 with $n=10$ , $max\mathrm{\_}try=100$ , and $max\mathrm{\_}itr=50$ . The curve tells us that in the first cycle of $j$ -loop, the inner for loop of Algorithm 1, no solution is found after $max\mathrm{\_}itr=50$ iterations of update of continuous assignment vector $\mathbf{s}$ . Then perturbation is given to $\mathbf{s}$ which causes a small jump of $L_{SU+\hat{\mathbf{c}}}$ at $itr=51$ and the second cycle of $j$ -loop starts and this time a solution is found after dozens of updates by thresholding $\mathbf{s}$ to a binary vector ${\mathbf{s}}_{*}$ .

The graph (c) in Figure 2 shows the scalability of computation time to find a solution up to $n=$ 10,000. We set $max\mathrm{\_}try=100$ , $max\mathrm{\_}itr=2000$ and plot the average of 10 measurements of time to find a solution. The graph seems to indicate good linearity w.r.t. $n$ up to $n=$ 10,000.

5.2. The HC Problem, Precomputation, and Another Solution Constraint

A HC is a cycle in a graph that visits every vertex exactly once and the HC problem is to determine if an HC exists in a given graph. It is an NP-complete problem and has been used as a programming example since the early days of ASP. Initially, it is encoded by a non-tight program containing positive recursion (Niemelä, 1999). Later a way of encoding by a program that is not tight but tight on its completion models is proposed (Lin & Zhao, 2003). We here introduce yet another encoding by a tight ground program inspired by SAT encoding proposed by Zhou (2020), where Zhou showed that the problem is solvable by translating six conditions listed in Figure 3 into a SAT problem. ¹²

Figure 3.

Conditions for Boolean satisfiability problem (SAT) encoding of a Hamiltonian cycle problem.

In what follows, we assume vertices are numbered from $1$ to $N$ = the number of vertices in a graph. We use $i\to j$ to denote an edge from vertex $i$ to vertex $j$ and $H_{i,j}$ to indicate there exists an edge from $i$ to $j$ in an HC. $U_{j,q}$ means vertex $j$ is visited at time $q$ ( $1\le j,q\le N$ ) and $\mathrm{one}\mathrm{\_}\mathrm{of}(a_1,\dots ,a_k)$ means that one of $\{a_1,\dots ,a_k\}$ is true. We translate these conditions into a program $P_3=\{(1),(2),(3)\}$ and constraints $K_3=\{(4),(5),(6)\}$ . To be more precise, the first condition $(1)$ is translated into a tight program just like a program $P_1$ (23). The conditions $\{(2),(3)\}$ constitute a tight definite program. Constraints $K_2=\{(4),(5),(6)\}$ are encoded as a set of implications of the form $\leftarrow L_1\wedge \cdots \wedge L_k$ , where $L_1,\dots ,L_k$ are literals. A set of $U_{j,q}$ atoms contained in a stable model of $P_3$ satisfying $K_3$ gives an HC.

We apply the above encoding to a simple HC problem for a graph $G_2$ in Figure 4. ¹³ There are six vertices and six HCs. ¹⁴ To solve this HC problem, we matricize $P_3$ and $K_3$ . There are 36 $H_{i,j}$ atoms ( $1\le i,j\le 6$ ) and 36 $U_{j,q}$ atoms ( $1\le j,q\le 6$ ). So there are 72 atoms in total. $P_3=\{(1),(2),(3)\}$ contains 197 rules in these 72 atoms and we translate $P_3$ into a pair of matrices $({\mathbf{C}}_3,{\mathbf{D}}_3)$ , where ${\mathbf{D}}_3$ is a $72\times 197$ binary matrix for disjunctions ¹⁵ and ${\mathbf{C}}_3$ is a $197\times 144$ matrix for conjunctions (rule bodies). Likewise $K_3=\{(4),(5),(6)\}$ is translated into a constraint matrix ${\hat{\mathbf{C}}}_{K_3}$ which is a $67\times 144$ binary matrix. A more detailed description of the encoding is available in the appendix (Appendix C). Then our task is to find a root $\mathbf{s}$ of $L_{SU+\hat{\mathbf{c}}}$ (10) constructed from these ${\mathbf{C}}_3$ , ${\mathbf{D}}_3$ , and ${\hat{\mathbf{C}}}_{K_3}$ in a 72 dimensional vector space by minimizing $L_{SU+\hat{\mathbf{c}}}$ to zero.

Figure 4.

A Hamiltonian cycle (HC) problem: (a) graph $G_2$ and (b) time and the number of different solutions.

We apply precomputation in the previous section to ${\mathbf{P}}_3=({\mathbf{C}}_3,{\mathbf{D}}_3)$ and ${\hat{\mathbf{C}}}_{K_3}$ to reduce program size. It takes 2.3 ms and detects 32 false stable atoms. It outputs a precomputed program ${\mathbf{P}}_3^{\prime }=({\mathbf{C}}_3^{\prime },{\mathbf{D}}_3^{\prime })$ and a constraint matrix ${\hat{\mathbf{C}}}_{K_3}^{\prime }$ of size ${\mathbf{D}}_3^{\prime }(40\times 90)$ , ${\mathbf{C}}_3^{\prime }(90\times 80)$ , and ${\hat{\mathbf{C}}}_{K_3}^{\prime }(52\times 80)$ , respectively, which is $1/4$ or $1/2$ of the original size. So precomputation removes $45\mathrm{\%}$ of atoms from the search space and returns much smaller matrices.

We run Algorithm 1 on ${\mathbf{P}}_3=({\mathbf{C}}_3,{\mathbf{D}}_3)$ with ${\hat{\mathbf{C}}}_{K_3}$ (no precomputation) and also on ${\mathbf{P}}_3^{\prime }=({\mathbf{C}}_3^{\prime },{\mathbf{D}}_3^{\prime })$ with ${\hat{\mathbf{C}}}_{K_3}^{\prime }$ (precomputation) using $max\mathrm{\_}try=20$ , $max\mathrm{\_}itr=200$ , and ${\ell }_2={\ell }_3=0.1$ and measure time to find a solution, that is, stable model satisfying constraints. The result is shown in Table (b) in Figure 4 as time(s), where time(s) is an average of 10 trials. The figures in the table, 2.08 s versus 0.66 s, ¹⁶ clearly demonstrate the usefulness of precomputation.

In addition to computation time, we examine the search power of different solutions in our approach by measuring the number of obtainable solutions. More concretely, we run Algorithm 1 seven times, and each time a stable model is obtained as a conjunction $L_1\wedge \cdots \wedge L_{72}$ of literals, we add a new constraint $\leftarrow L_1\wedge \cdots \wedge L_{72}$ to previous constraints, thereby forcibly computing a new stable model in next trial. We call such use of constraint another solution constraint. Since there are at most six solutions, the number of solutions obtained by seven trials is at most six. We repeat a batch of seven trials 10 times and take the average of the number of solutions obtained by each batch. The average is denoted as #solutions in Table (b) which indicates that 5.7 solutions out of 6, almost all solutions, are obtained by seven trials using another solution constraint.

Summing up, the figures in Table (b) exemplify the effectiveness of precomputation which significantly reduces computation time and returns a more variety of solutions when combined with another solution constraint.

5.3. LF Heuristics and Precomputation on Loopy Programs

So far we have been dealing with tight programs which have no loop and hence have no loop formulas. We here deal with non-tight programs containing loops and examine how LF heuristics, $L{F}_{max}$ and $L{F}_{min}$ , introduced in the previous section work. We use an artificial non-tight program $P_{4\mathrm{\_}n}$ (with no constraint) shown below that has exponentially many loops.

$$\begin{array}{rl}{P}_{4\mathrm{\_}n}& =\{\begin{array}{lll}{a}_0& \leftarrow & a_1\wedge \cdots \wedge a_n\\ a_0& \leftarrow & \mathrm{\neg }{a}_{n+1}\\ \dots \\ a_{2i-1}& \leftarrow & a_0\vee a_{2i}& \text{for }i:1\le i\le n/2\\ a_{2i}& \leftarrow & a_0\vee a_{2i-1}& \text{for }i:1\le i\le n/2\\ \dots \\ a_{n+1}& \leftarrow & a_{n+1}\end{array}\end{array}$$

We here consider an even $n$ , then $P_{4\mathrm{\_}n}$ program has $n+2$ atoms $\{a_0,a_1,\dots ,a_n,a_{n+1}\}$ , $2^{n/2}+1$ supported models and one stable model $\{a_0,a_1,\dots ,a_n\}$ . There are $n/2+1$ minimal loops $\{a_1,a_2\},\dots ,\{a_{n-1},a_n\},\{a_{n+1}\}$ and a maximal loop $\{a_0,a_1,\dots ,a_n\}$ . The set of loop formulas for LF heuristics are computed as follows:

$$\begin{array}{rl}L{F}_{max}& =\{(a_0\wedge a_1\wedge \cdots \wedge a_n)\to \mathrm{\neg }{a}_{n+1},a_{n+1}\to \perp \},\\ L{F}_{min}& =\{(a_1\wedge a_2)\to a_0,\dots ,(a_{n-1}\wedge a_n)\to a_0,a_{n+1}\to \perp \}.\end{array}$$

Note that although there are $2^{n/2}+1$ supported models, there is only one stable model. So $L{F}_{max}$ and $L{F}_{min}$ are expected to exclude $2^{n/2}$ supported models.

After translating $P_{4\mathrm{\_}n}$ into a matricized program ${\mathbf{P}}_{4\mathrm{\_}n}=({\mathbf{C}}_{4\mathrm{\_}n},{\mathbf{D}}_{4\mathrm{\_}n})$ , where ${\mathbf{C}}_{4\mathrm{\_}n}$ is a $(2n+3)\times (2n+4)$ binary matrix and ${\mathbf{D}}_{4\mathrm{\_}n}$ is a $(n+2)\times (2n+3)$ binary matrix, respectively, we compute a stable model of $P_{4\mathrm{\_}n}$ for various $n$ by Algorithm 1 that minimizes $L_{SU+\hat{\mathbf{c}}+LF}$ (19) with coefficient ${\ell }_3=0$ for the constraint term (because of no use of constraints) using Jacobian $J_{a_{SU+\hat{\mathbf{c}}+LF}}$ (11).

Below is an example of the program $P_{4\mathrm{\_}n}$ , where $n=4$ .

$$\begin{array}{rl}{P}_{4\mathrm{\_}4}& =\begin{array}{l}\{\begin{array}{lll}{a}_0& \leftarrow & a_1\wedge a_2\wedge a_3\wedge a_4\\ a_0& \leftarrow & \mathrm{\neg }{a}_5\\ a_1& \leftarrow & a_0& a_1\leftarrow a_2\\ a_2& \leftarrow & a_0& a_2\leftarrow a_1\\ a_3& \leftarrow & a_0& a_3\leftarrow a_4\\ a_4& \leftarrow & a_0& a_4\leftarrow a_3\\ a_5& \leftarrow & a_5\end{array}\end{array}\end{array}$$

This program has three minimal loops $\{a_1,a_2\},\{a_3,a_4\},a_5$ , and a maximal loop $\{a_0,a_1,a_2,a_3,a_4\}$ . There are 11 rules and six atoms, so ${\mathbf{C}}_{4\mathrm{\_}4}$ is a $(11\times 12)$ binary matrix.

Since all supported models of $P_{4\mathrm{\_}n}$ except for one stable model are non-stable, even if $L{F}_{max}$ and $L{F}_{min}$ are used to guide the search process towards a stable model, Algorithm 1 is likely to return a non-stable model. We can avoid such a situation by the use of another solution constraint.

To verify it, we examine the pure effect of another solution constraint that guides the search process to compute a model different from previous ones. Without using $L{F}_{max}$ or $L{F}_{min}$ heuristics, we repeatedly run Algorithm 1 with/without another solution constraint for ${10}^4$ trials with $n=4$ , $max\mathrm{\_}try=20$ , $max\mathrm{\_}itr=50$ , ${\ell }_2={\ell }_3=0.1$ and measure time to find a stable model and count the number of trials until then. We repeat this experiment 10 times and take the average. The result is shown in Table 2.

Table 2.

The Effect of Another Solution Constraint.

Another solution constraint	Time(s)	# trials
Not used	11.46 (0.41)	${10}^4(0)$
Used	0.09 (0.13)	3.5 (1.6)

The figure ${10}^4(0)$ in Table 2 in the case of no use of another solution constraint means Algorithm 1 always exhausts ${10}^4$ trials without finding a stable model (due to implicit bias in Algorithm 1). When another solution constraint is used, however, it finds a stable model in 0.09 s after 3.5 trials on average. Thus Table 2 demonstrates the necessity and effectiveness of another solution constraint to efficiently explore the search space.

We next compare the effectiveness of LF heuristics and that of precomputation under another solution constraint. For $n=10,\dots ,50$ , we repeatedly run Algorithm 1 using $L_{SU+\hat{\mathbf{c}}+LF}$ with $max\mathrm{\_}try=10,max\mathrm{\_}itr=100$ on matricized ${\mathbf{P}}_{4\mathrm{\_}n}=({\mathbf{C}}_{4\mathrm{\_}n},{\mathbf{D}}_{4\mathrm{\_}n})$ (and no constraint matrix) to compute supported (stable) models. Coefficients in $L_{SU+\hat{\mathbf{c}}+LF}$ are set to ${\ell }_2=0.1,{\ell }_3=0,\mathrm{and}{\ell }_4=1$ . To be more precise, for each $n$ and each case of $L{F}_{max}$ , $L{F}_{min}$ , precomputation (without $\{L{F}_{max},L{F}_{min}\}$ ) and no $\{L{F}_{max},L{F}_{min},\text{precomputation}\}$ , we run Algorithm 1 at most 100 trials to measure time to find a stable model and count the number of supported models computed till then. We repeat this computation 10 times and take the average and obtain graphs in Figure 5.

Figure 5.

The effect of loop formula (LF) heuristics and precomputation on program $P_{4\mathrm{\_}n}$ : (a) time to find a stable model and (b) #computed models.

In Figure 5, no_LF means no use of $\{L{F}_{max},L{F}_{min}\}$ heuristics. Also no_LF_pre means no_LF is applied to precomputed $P_{4\mathrm{\_}n}$ . ¹⁷

We can see from graph (a) in Figure 5 that computation time is $L{F}_{min}>L{F}_{max}>\text{no\_LF}>\text{no\_LF\_pre}$ . This means that using LF heuristics is not necessarily a good policy. They might cause extra computation to reach the same model. Concerning the number of non-stable models computed redundantly, graph (b) in Figure 5 tells us that $L{F}_{min}$ allows computing redundant non-stable models but the rest, $L{F}_{max}$ , no_LF, and no_LF_pre, return a stable model without computing redundant non-stable models. This shows first that $L{F}_{max}$ works correctly to suppress the computation of non-stable models and second that the $L{F}_{min}$ heuristics works adversely, that is, guiding the search process away from the stable model. This somewhat unexpected result indicates the need of (empirical) choice of LF heuristics.

Finally, to examine the effectiveness of precomputation more precisely, we apply precomputation to a more complex program $P_{5\mathrm{\_}nk}$ . It is a modification of $P_{4\mathrm{\_}n}$ by adding self-loops of $k$ atoms as illustrated by (a) in Figure 6. The addition of self-loop causes the choice of $a_{n+j}$ ( $1\le j\le k$ ) being true or being false in the search process. $P_{5\mathrm{\_}nk}$ has $(2^{n/2}-1)(2^k-1)+1$ supported models but has just one stable model $\{a_0,a_1,\dots ,a_n\}$ .

Figure 6.

Precomputation applied to program $P_{5\mathrm{\_}nk}$ : (a) a non-tight program $P_{5\mathrm{\_}nk}$ and (b) scalability of precomputation w.r.t. $P_{5\mathrm{\_}nk}$ .

We compute a stable model by running Algorithm 1 on precomputed $P_{5\mathrm{\_}nk}$ without using LF heuristics up to $n=k=5000$ . When precomputation is applied to $P_{5\mathrm{\_}nk}$ , where $n=k=5000$ , it detects $5000$ false stable atoms and downsizes the matrices in ${\mathbf{P}}_{5\mathrm{\_}nk}=({\mathbf{C}}_{5\mathrm{\_}nk},{\mathbf{D}}_{5\mathrm{\_}nk})$ from ${\mathbf{D}}_{5\mathrm{\_}nk}(10001\times 15002)$ to ${\mathbf{D}}_{5\mathrm{\_}nk}^{\prime }(5001\times 10002)$ and from ${\mathbf{C}}_{5\mathrm{\_}nk}(15002\times 20002)$ to ${\mathbf{C}}_{5\mathrm{\_}nk}^{\prime }(10002\times 10002)$ . Thus precomputed ${\mathbf{P}}_{5\mathrm{\_}nk}^{\prime }=({\mathbf{C}}_{5\mathrm{\_}nk}^{\prime },{\mathbf{D}}_{5\mathrm{\_}nk}^{\prime })$ is downsized to 1/3 of the original $P_{5\mathrm{\_}nk}$ .

We run Algorithm 1 on $P_{5\mathrm{\_}nk}^{\prime }$ with ${\ell }_2={\ell }_3=0.1$ and $max\mathrm{\_}try=10$ , $max\mathrm{\_}itr=100$ at most 100 trials to measure time to find a stable model ten times for each $n=1000,\dots ,5000$ and take the average. At the same time, we also run clingo (version 5.6.2) on $P_{5\mathrm{\_}nk}$ and similarly measure time. Graph (b) in Figure 6 is the result. It shows that as far as computing a stable model of $P_{5\mathrm{\_}nk}$ is concerned, our approach comes close to clingo. However, this is due to a very specific situation that precomputation removes all false atoms $\{a_{n+1},\dots ,a_{n+k}\}$ in the stable model of $P_{5\mathrm{\_}nk}$ and Algorithm 1 run on the precomputed ${\mathbf{P}}_{5\mathrm{\_}nk}^{\prime }=({\mathbf{C}}_{5\mathrm{\_}nk}^{\prime },{\mathbf{D}}_{5\mathrm{\_}nk}^{\prime })$ detects the stable model only by thresholding $\mathbf{s}$ before starting any update of $\mathbf{s}$ . So what graph (b) really suggests seems to be the importance of optimization of a program like precomputation, which is to be developed further in our approach.

6. Related Work

The most closely related work is Aspis et al. (2020) and Takemura and Inoue (2022). As mentioned in Section 1, our approach differs from them in three points: (1) theoretically, the exclusion of non-stable models by loop formulas, (2) syntactically, no restriction on acceptable programs, and (3) practically, incorporation of constraints. Concerning performance, they happen to use the same $N$ -negative loops program which consists of $N$ copies (alphabetic variants) of a program $\{a\leftarrow \mathrm{\neg }b,b\leftarrow \mathrm{\neg }a\}$ . According to Aspis et al. (2020), the success rate w.r.t. $N$ of returning a supported model goes from one initially to almost zero at $N=64$ by Aspis et al. (2020), while it keeps one till $N=100$ by Takemura and Inoue (2022). We tested the same program with $max\mathrm{\_}try=20$ , $max\mathrm{\_}itr=100$ , and observed that the success rate keeps one till $N=$ 10,000.

Although our approach is non-probabilistic, that is, purely linear algebraic, there are probabilistic differentiable approaches for ASP. Differentiable ASP/SAT (Nickles, 2018) iteratively samples a stable model by an ASP solver a la ASSAT (Lin & Zhao, 2004). The solver decides the next decision literal based on the derivatives of a cost function which is the MSE between the target probabilities and predicted probabilities computed from the sampled stable models via parameters associated with “parameter atoms” in a program.

NeurASP (Yang et al., 2020) uses an ASP solver to obtain stable models including “neural atoms” for a program. They are associated with probabilities learned by deep learning and the likelihood of an observation (a set of ASP constraints) is computed from them. The whole learning is carried out by backpropagating the likelihood to neural atoms to parameters in a neural network.

Similarly to NeurASP, SLASH (Skryagin et al., 2023, 2022) uses an ASP solver to compute stable models for a program containing “neural probabilistic predicates.” Their probabilities are dealt with by neural networks and probabilistic circuits. The latter makes it possible to compute a joint distribution of the class category and data. Both NeurASP and SLASH are examples of symbolic ASP solver-based neuro-symbolic systems, where they include a neural frontend to process the perception part of the problem, and a symbolic backend which typically is the ASP solver. Therefore, the neural frontend does not need to be involved in the computational details and problems associated with computing stable models (Section 4).

Independently of ASP solver-based approaches mentioned above, Sato and Kojima proposed a differentiable approach to sampling supported models of (non-propositional) probabilistic normal logic programs (Sato & Kojima, 2019, 2020). They encode programs by matrices and formulate the problem of sampling supported models as repeatedly computing a fixedpoint of some differentiable equations. They solve the equations in vector spaces by minimizing a non-negative cost function defined by Frobenius norm. More recently, Takemura and Inoue (2024) proposed a neuro-symbolic learning pipeline for distant supervision tasks, which leverages differentiable computation of supported models. Similarly to this work, they encode normal logic programs into matrices and define a differentiable loss function which is based on the supported model semantics.

As for the non-differentiable linear-algebraic approaches to logic programming, Nguyen et al. adopted matrix encoding for propositional normal logic programs based on Sakama et al. (2017) and proposed to compute stable models in vector spaces by a generate-and-test approach using sparse representation (Nguyen et al., 2022).

6.1. Connection to Neural Network Computation

At this point, it is quite interesting to see the connection of our approach to neural network computation. In (6), we compute $\mathbf{M}$ and $\mathbf{d}=\mathbf{D}\mathbf{M}$ . We point out that the computation of this $\mathbf{d}$ is nothing but the output of a forward process of a single layer ReLU network from an input vector $\mathbf{s}$ . Consider the computation of $\mathbf{M}=(\mathbf{1}-{min}_1(\mathbf{C}(\mathbf{1}-{\mathbf{s}}^{\delta })))$ . We rewrite this using $1-\min (x,1)=\mathrm{ReLU}(1-x)$ to

$$\begin{array}{rl}\mathbf{M}& =\mathbf{1}-{min}_1({\mathbf{C}}^{pos}(\mathbf{1}-\mathbf{s})+{\mathbf{C}}^{neg}\mathbf{s})\\ & =\mathrm{ReLU}(\mathbf{W}\mathbf{s}+\mathbf{b}),\\ & \text{where}\mathbf{C}=[{\mathbf{C}}^{pos}{\mathbf{C}}^{neg}],\mathbf{W}={\mathbf{C}}^{pos}-{\mathbf{C}}^{neg},\mathbf{b}=\mathbf{1}-{\mathbf{C}}^{pos}\mathbf{1}.\end{array}$$

So $\mathbf{M}$ is the output of a ReLU network having a weight matrix $\mathbf{W}={\mathbf{C}}^{pos}-{\mathbf{C}}^{neg}$ and a bias vector $\mathbf{b}=\mathbf{1}-{\mathbf{C}}^{pos}\mathbf{1}$ . Then ${\min }_1(\mathbf{d})={\min }_1(\mathbf{D}\mathbf{M})={\min }_1(\mathbf{D}\cdot \mathrm{ReLU}(\mathbf{W}\mathbf{s}+\mathbf{b}))$ is the output of a ReLU network with a single hidden layer and a linear output layer represented by $\mathbf{D}$ having ${\min }_1(\cdot )$ as activation function. Also when we compute a supported model $\mathbf{s}$ , we minimize $L_{SU+\hat{\mathbf{c}}}$ (6) which contains an MSE error term $L_{sq}=\Vert {\min }_1(\mathbf{d})-\mathbf{s}{\Vert }^2$ using $J_{a_{SU+\hat{\mathbf{c}}}}$ (11). This is precisely back propagation from learning data $\mathbf{s}$ .

Thus we may say that our approach is an integration of ASP semantics and neural computation and provides a neuro-symbolic (Sarker et al., 2021) way of ASP computation. Nonetheless, there is a big difference. In standard neural network architecture, a weight matrix $\mathbf{W}$ and a bias vector $\mathbf{b}$ are independent. In our setting, they are interdependent and they faithfully reflect the logical structure of a program.

7. Conclusion

We proposed an end-to-end approach for computing stable models satisfying given constraints. We matricized a program and constraints and formulated stable model computation as a minimization problem in vector spaces of a non-negative cost function. We obtain a stable model satisfying constraints as a root the cost function by Newton’s method.

By incorporating all loop formula constraints introduced in Lin-Zhao’s theorem (Lin & Zhao, 2004) into the cost function to be minimized, we can prevent redundant computation of non-stable models, at the cost of processing exponentially many loop formulas. Hence, we introduced precomputation which downsizes a program while preserving stable model semantics and also two heuristics that selectively use loop formulas. Then we confirmed the effectiveness of our approach including precomputation and loop formula heuristics by simple examples.

Future work could focus on improving the integration of neural networks with this proposed end-to-end approach to tackle neuro-symbolic benchmark tasks that require both perception and reasoning. We also aim to improve the optimization techniques, such as precomputation, to enhance efficiency and scalability.