¹¹institutetext: National Institute of Informatics, Japan
¹¹email: [email protected] ²²institutetext: National Institute of Informatics, Japan
²²email: [email protected] ³³institutetext: National Institute of Informatics, Japan
³³email: [email protected]

Towards end-to-end ASP computation

Taisuke Sato 11 0000-0001-9062-0729 Akihiro Takemura 22 0000-0003-4130-8311 Katsumi Inoue 33 0000-0002-2717-9122

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 [14] 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 3-coloring and Hamiltonian cycle problems. As our approach is purely numerical and only contains vector/matrix operations, acceleration by parallel technologies such as many-cores and GPUs is expected.

Keywords:

Answer Set Programming end-to-end ASP vector space cost minimization loop formula supported model stable model

1 Introduction

Computing stable model semantics [9] lies at the heart of answer set programming (ASP) [19, 16, 11] and there have been a variety of approaches proposed so far. Early approaches such as Smodels [20] used backtracking. Then the concept of loop formula was introduced and approaches that use a SAT solver to compute stable models based on Lin-Zhao’s theorem [14] were proposed. They include ASSAT [14] and Cmodels [10] for example. Later more elaborated approaches such as Clasp [7, 8] 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 [18] computes stable models by an ASP solver that utilizes derivatives of a cost function. More recently NeurASP [29] and SLASH [26] 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. A step towards end-to-end computation was taken by Aspis et al. [2] and Takemura and Inoue [27]. They formulated the computation of supported models, a super class of stable models¹¹1 Supported models and stable models of a propositional normal logic program coincide when the program is tight (no infinite call chain through positive goals) [5, 4]. , entirely as fixedpoint computation in vector spaces and obtain supported models represented by binary vectors [1, 17]. Also in the context of probabilistic modeling, an end-to-end sampling of supported models was proposed by Sato and Kojima in [24]. However there still remains a gap between computing stable models and computing supported ones.

In this paper, 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 [14] together with constraints directly in vector spaces as a cost minimization problem, thereby no use of symbolic ASP or SAT solvers involved. Also as our approach is purely numerical and only contains vector/matrix operations, acceleration by parallel technologies such as many-cores and GPUs is expected.

Technically, Lin-Zhao’s theorem [14] states that a stable model of a ground normal logic program coincides with a supported model which satisfies “loop formulas” associated with the program. They 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 [12]. 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. In a symbolic approach, the latter would mean allowing non-stable model computation but in our continuous approach, it means gracefully degraded gradient information for the continuous search process.

Our end-to-end computing framework differs from those by [2] and [27] in that they basically compute supported models and the computing process itself has no mechanism such as loop formulas to exclude non-stable models. Also we impose no restriction on the syntax of programs like the MD condition [2] or the SD condition [27], any propositional normal logic program is acceptable. More importantly we incorporate the use of constraints, i.e. rules like $\leftarrow a\wedge\neg b$ , which make ASP programming smooth and practical.

Hence our contributions include

•

a proposal of end-to-end computing of stable models in vector spaces for unrestricted 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 paper is to establish theoretical feasibility of end-to-end ASP computing in vector spaces, programming examples are small and implementation is of preliminary nature.

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 3-coloring problem and the HC problem. We there compare performance of precomputation and loop formula heuristics. Section 6 contains related work and Section 7 is conclusion.

2 Preliminaries

In this paper, 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²²2 We equate propositional variables with atoms. called the head, $G$ a conjunction of literals³³3 A literal is an atom (positive literal) or its negation (negative literal). called the body of the rule respectively. We suppose $P$ is written in a given set of atoms ${\cal A}$ but usually assume ${\cal A}$ = atom( $P$ ), i.e. the set of atoms occurring in $P$ . We use $G^{+}$ to denote the conjunction of positive literals in $G$ . $G$ may be empty. The empty conjunction is always true. We call $a\leftarrow G$ rule for $a$ . Let $a\leftarrow G_{1},\ldots,a\leftarrow G_{m}$ be rules for $a$ in $P$ . When $m>0$ , put $\mbox{iff}(a)=a\Leftrightarrow G_{1}\vee\cdots\vee G_{m}$ . When $m=0$ , i.e. there is no rule for $a$ , put $\mbox{iff}(a)=a\Leftrightarrow\bot$ where $\bot$ is a special symbol representing the empty disjunction which is always false. We call $\mbox{iff}(a)$ the completed rule for $a$ . The completion of $P$ , comp( $P$ ), is defined as comp( $P$ ) = $\{\mbox{iff}(a)\mid\mbox{atom $a$ 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$ .

A model (assignment) $I$ over a set of atoms ${\cal A}$ is a mapping which determines the truth value of each atom $a\in{\cal 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},\ldots,F_{k}\}$ by considering $F$ as a conjunction $F_{1}\wedge\cdots\wedge F_{k}$ . For convenience, we always equate $I$ with $\{a\in{\cal A}\mid I\models a\}$ , i.e. the set of atoms true in $I$ . When $I$ satisfies all rules in the program $P$ , i.e. $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{\cal A}\mid P\vdash a\}$ , i.e. the set of atoms provable from $P$ . A model $I$ of comp( $P$ ) is called a supported model of $P$ [17]. When $P$ is a definite program, its least model is also a supported model. In general, for non-definite $P$ , there can be zero or multiple supported models. 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) reduction 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$ [9]. $P$ may have zero or multiple stable models like supported models. Since the existence of a stable model is NP-complete [16] and so is a supported model, their computation is expected to be hard.

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\leq i\leq h$ ) is a conjunction of literals and called a disjunct of $F$ . When $F$ has no disjunct, $F$ is false. $F$ is called full when every $d_{i}$ is a conjunction of $n$ distinct literals.

A walk in a directed graph is a sequence $v_{1}\rightarrow v_{2}\rightarrow\cdots\rightarrow v_{k}$ ( $k\geq 1$ ) of vertices representing the corresponding non-zero sequence of edges $(v_{1},v_{2}),\ldots,(v_{k-1},v_{k})$ . When $v_{k}=v_{1}$ , it is said to be closed. A cycle is a closed walk $v_{1}\rightarrow v_{2}\rightarrow\cdots\rightarrow v_{k}\rightarrow v_{1}$ where $\{v_{1},\ldots,v_{k}\}$ are all distinct. A Hamiltonian cycle is a cycle which visits every vertex 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 $\mbox{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 [5]⁴⁴4 In [5], it is called “positive-order-consistent”. when $\mbox{pdg}(P)$ is acyclic, i.e. has no cycle. A loop $L=\{a_{1},\ldots,a_{k}\}$ $(k>0)$ in $P$ is a set of atoms where for any pair of atoms $a_{1}$ and $a_{2}$ in $L$ ( $a_{1}=a_{2}$ allowed), there is a path in $\mbox{pdg}(P)$ from $a_{1}$ to $a_{2}$ and also from $a_{2}$ to $a_{1}$ .

We denote vectors by bold italic 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. $\|{\mathbf{a}}\|_{1}$ and $\|{\mathbf{a}}\|_{2}$ respectively denote the 1-norm and 2-norm of ${\mathbf{a}}$ where $\|{\mathbf{a}}\|_{1}=\sum|{\mathbf{a}}(i)|$ and $\|{\mathbf{a}}\|_{2}=\sqrt{\sum{\mathbf{a}}(i)^{2}}$ . We use ${\mathbf{1}}$ to denote an all-ones vector of appropriate dimension. A model $I$ over a set ${\cal A}=\{a_{1},\ldots,a_{n}\}$ of $n$ ordered atoms is equated with an $n$ -dimensional binary vector ${\mathbf{u}}_{I}\in\{1,0\}^{n}$ such that ${\mathbf{u}}_{I}(i)=1$ if $a_{i}$ is true in $I$ and ${\mathbf{u}}_{I}(i)=0$ otherwise ( $1\leq i\leq n$ ). ${\mathbf{u}}_{I}$ is called the vectorized $I$ .

Italic capital letters such as ${A}$ stand for a matrix. We use ${A}(i,j)$ to denote the $i,j$ -th element of ${A}$ , ${A}(i,:)$ the $i$ -th row of ${A}$ and ${A}(:,j)$ the $j$ -th column of ${A}$ , respectively. We often consider one dimensional matrices as (row or column) vectors. $\|{A}\|_{F}$ denotes the Frobenius norm of ${A}$ . Let ${A},{B}\in\mathbb{R}^{m\times n}$ be $m\times n$ matrices. ${A}\odot{B}$ denotes their Hadamard product, i.e., $({A}\odot{B})(i,j)={A}(i,j){B}(i,j)$ for $i,j(1\leq i\leq m,1\leq j\leq n)$ . $[{A};{B}]$ designates the $2m\times n$ matrix of ${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 $\mbox{min}_{1}(x)=\mbox{min}(x,1)$ that returns the lesser of 1 and $x$ as an activation function which is related to the popular activation function $\mbox{ReLU}(x)=\mbox{max}(x,0)$ by $1-\mbox{min}_{1}(x)=\mbox{ReLU}(1-x)$ . $\mbox{min}_{1}({A})$ denotes the result of component-wise application of $\mbox{min}_{1}(x)$ to matrix ${A}$ . We also introduce thresholding notation. Suppose $\theta$ is a real number and ${\mathbf{a}}$ an $n$ -dimensional vector. Then $[{\mathbf{a}}\leq\theta]$ denotes a binary vector obtained by thresholding ${\mathbf{a}}$ at $\theta$ where for $i(1\leq i\leq n)$ , $[{\mathbf{a}}\leq\theta](i)=1$ if ${\mathbf{a}}(i)\leq\theta$ and $[{\mathbf{a}}\leq\theta](i)=0$ otherwise. $[{\mathbf{a}}\geq\theta]$ is treated similarly. We extend thresholding to matrices. Thus $[A\leq 1]$ means a matrix such that $[A\leq 1](i,j)=1$ if $A(i,j)\leq 1$ and $[A\leq 1](i,j)=0$ otherwise. For convenience, we generalize bit inversion to an $n$ -dimensional vector ${\mathbf{a}}$ and use an expression $1-{\mathbf{a}}$ to denote the $n$ -dimensional vector such that $(1-{\mathbf{a}})(i)=1-{\mathbf{a}}(i)$ for $i\,(1\leq i\leq n)$ . $1-{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

First we encode programs as matrices. For concreteness, we explain by an example (generalization is easy). Suppose we are given a program $P_{0}$ below containing three rules $\{r_{1},r_{2},r_{3}\}$ in a set of atoms ${\cal A}_{0}=\{p,q,r\}$ .

\displaystyle P_{0}

\displaystyle=

\displaystyle\begin{array}[]{ll}\left\{\;\begin{array}[]{lllll}p&\leftarrow&q\;\wedge\,\neg r&&\mbox{: rule $\,r_{1}\,$ for $p$}\\ p&\leftarrow&\neg q&&\mbox{: rule $\,r_{2}\,$ for $p$}\\ q&&&&\mbox{: rule $\,r_{3}\,$ for $q$}\\ \end{array}\right.\end{array}

(5)

Assuming atoms are ordered as $p,q,r$ and correspondingly so are the rules $\{r_{1},r_{2},r_{3}\}$ as in (5), we encode $P_{0}$ as a pair of matrices $(D_{0},Q_{0})$ . Here $Q_{0}$ represents conjunctions (the bodies of $\{r_{1},r_{2},r_{3}\}$ ) and $D_{0}$ their disjunctions so that they jointly represent $P_{0}$ .

\displaystyle Q_{0}\;

\displaystyle=

\displaystyle\begin{array}[]{ll}\begin{matrix}\;\;\;\;p&\;q&\;r&\!\!\!\!\neg p&\!\!\!\neg q&\!\!\!\!\neg r\\ \vspace{-0.8em}\end{matrix}\\ \begin{bmatrix}\;\;0&\,1&\,0&\;0&\,0&\,1\;\;\\ \;\;0&\,0&\,0&\;0&\,1&\,0\;\;\\ \;\;0&\,0&\,0&\;0&\,0&\,0\;\;\\ \end{bmatrix}&\begin{matrix}\;\;\mbox{: $r_{1}$ has the body $q\wedge\neg r$}\hfill\\ \;\;\mbox{: $r_{2}$ has the body $\neg q$}\hfill\\ \;\;\mbox{: $r_{3}$ has the empty body}\hfill\\ \end{matrix}\end{array}

(8)

\displaystyle D_{0}\;

\displaystyle=

\displaystyle\begin{array}[]{ll}\begin{matrix}\;\;\;\;r_{1}&r_{2}&\!r_{3}\\ \vspace{-0.8em}\end{matrix}\\ \begin{bmatrix}\;\;1&\,1&\,0\;\;\\ \;\;0&\,0&\,1\;\;\\ \;\;0&\,0&\,0\;\;\\ \end{bmatrix}&\begin{matrix}\;\;\mbox{: $p\,$ has two rules $r_{1}$ and $r_{2}$}\hfill\\ \;\;\mbox{: $q\;$ has one rule $r_{3}$}\hfill\\ \;\;\mbox{: $r\;$ has no rule}\hfill\\ \end{matrix}\end{array}

(11)

As can be seen, $Q_{0}$ represents conjunctions in $P_{0}$ in such a way that $Q_{0}(1,:)$ for example represents the conjunction $q\wedge\neg r$ of the first rule in $P_{0}$ by setting $Q_{0}(1,2)=Q_{0}(1,6)=1$ and so on. $D_{0}$ represents disjunctions of rule bodies. So $D_{0}(1,1)=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\neg r)\vee\neg q$ for $p$ .

In general, a program $P$ that has $m$ rules in $n$ atoms is numerically encoded as a pair $P=(D,Q)$ of binary matrices $D\in\{0,1\}^{n\times m}$ and $Q\in\{0,1\}^{m\times 2n}$ , which we call a matricized $P$ . $Q$ represents rule bodies in $P$ . Suppose atoms are ordered like ${\cal A}=\{a_{1},\ldots,a_{n}\}$ and similarly rules are ordered like $\{r_{1}:a_{i_{1}}\leftarrow G_{1},\ldots,r_{m}:a_{i_{m}}\leftarrow G_{m}\}$ . Then the $j$ -th row $Q(j,:)$ ( $1\leq j\leq 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\neg a_{i_{p+1}}\wedge\cdots\wedge\neg a_{i_{p+q}}$ ( $1\leq p,q\leq n$ ). Then $Q(j,:)$ is a zero vector except for $Q(j,i_{1})=\cdots=Q(j,i_{p})=Q(j,n+i_{p+1})=\cdots=Q(j,n+i_{p+q})=1$ . $D$ combines these conjunctions as a disjucntion (DNF) for each atom in ${\cal A}$ . If the $i$ -th atom $a_{i}\in{\cal A}$ ( $1\leq i\leq n$ ) has rules $\{a_{i}\leftarrow G_{j_{1}},\ldots,a_{i}\leftarrow G_{j_{s}}\}$ in $P$ , we put $D(i,j_{1})=\cdots=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}$ : $\mbox{iff}(a_{i})=a_{i}\Leftrightarrow G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ . If $a_{i}$ has no rule, we put $D(i,j)=0$ for all $j$ ( $1\leq j\leq m$ ). Thus the matricized $P=(D,Q)$ can represent the completed program $\mbox{comp}(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 ${\cal A}$ of atoms. Recall that $I$ is equated with a subset of ${\cal A}$ . We inductively define the relation “a formula $F$ is true in $I$ ”⁵⁵5 Equivalently, $I$ satisfies $F$ . , $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\neg F$ iff $I\not\models F$ . When $F$ is a disjunction $F_{1}\vee\cdots\vee F_{k}$ ( $k\geq 0$ ), $I\models F$ iff there is some $i$ ( $1\leq i\leq k$ ) s.t. $I\models F_{i}$ . So the empty disjunction ( $k=0$ ) is always false. We consider a conjunction $F_{1}\wedge\cdots\wedge F_{k}$ as a syntax sugar for $\neg(\neg F_{1}\vee\cdots\vee\neg F_{k})$ 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\neg a_{i_{p+1}}\wedge\cdots\wedge\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}},\ldots,a_{i_{p}}\}\subseteq I$ and $\{a_{i_{p+1}},\ldots,a_{i_{p+q}}\}\cap I=\emptyset$ . Also let $\mbox{iff}(a_{i})=a_{i}\Leftrightarrow G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ be the completed rule for an atom $a_{i}$ in $P$ . We see $I\models\mbox{iff}(a)$ iff $\left(a_{i}\in I\;\mbox{iff}\;I\models G_{j_{1}}\vee\cdots\vee G_{j_{s}}\right)$ .

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

Consider a matricized program $P=(D,Q)$ $(D\in\{0,1\}^{n\times m}$ , $Q\in\{0,1\}^{m\times 2n})$ and its $j$ -th rule $r_{j}$ having a body $G_{j}$ represented by $Q(j,:)$ . Compute $Q(j,:){\mathbf{u}}_{I}^{\delta}$ which is the number of true literals in $I$ in $G_{j}$ and compare it with the number of literals $|Q(j,:)|_{1}$ ⁶⁶6 $|{\mathbf{v}}|_{1}=\sum_{i}|{\mathbf{v}}(i)|$ is the 1-norm of a vector ${\mathbf{v}}$ . in $G_{j}$ . When $|Q(j,:)|_{1}=Q(j,:){\mathbf{u}}_{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 $Q(j,:)(1-{\mathbf{u}}_{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 $\mbox{min}_{1}(Q(j,:)(1-{\mathbf{u}}_{I}^{\delta}))$ which is $1$ if there is a false literal, and $0$ otherwise. Consequently $1-\mbox{min}_{1}(Q(j,:)(1-{\mathbf{u}}_{I}^{\delta}))=1$ if there is no false literal in $G_{j}$ and vice versa. In other words, $1-\mbox{min}_{1}(Q(j,:)(1-{\mathbf{u}}_{I}^{\delta}))$ computes $I\models G_{j}$ .

Now let $\{a_{i}\leftarrow G_{j_{1}},\ldots,a_{i}\leftarrow G_{j_{s}}\}$ be an enumeration of rules for $a_{i}\in{\cal A}$ and $G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ the disjunction of the rule bodies. $d_{i}=\sum_{t=1}^{s}(1-\mbox{min}_{1}(Q(j_{t},:)(1-{\mathbf{u}}_{I}^{\delta})))$ is the number of rule bodies in $\{G_{j_{1}},\ldots,G_{j_{s}}\}$ that are true in $I$ . Noting $D(i,j)=1$ if $j\in\{j_{1},\ldots,j_{s}\}$ and $D(i,j)=0$ otherwise by construction of $D$ in $P=(D,Q)$ , we replace the summation $\sum_{t=1}^{s}$ by matrix multiplication and obtain $d_{i}=D(i,:)(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}_{I}^{\delta})))$ . Introduce a column vector ${\mathbf{d}}_{I}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}_{I}^{\delta})))$ . We have ${\mathbf{d}}_{I}(i)=d_{i}$ = the number of rules for $a_{i}$ whose bodies are true in $I$ ( $1\leq i\leq n$ ).

Proposition 1

Let $P=(D,Q)$ be a matricized program $P$ in a set of atoms ${\cal A}$ and ${\mathbf{u}}_{I}$ a vectorized model $I$ over ${\cal A}$ . Put ${\mathbf{d}}_{I}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}_{I}^{\delta})))$ . It holds that

\displaystyle I\models\mbox{comp}(P)

if-and-only-if

\displaystyle\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})\|_{2}=0.

(12)

(Proof) Put $n=|{\cal A}|$ . Suppose $I\models\mbox{comp}(P)$ and write $\mbox{iff}(a_{i})$ , the completed rule for an atom $a_{i}\in{\cal A}$ ( $1\leq i\leq n$ ), as $\mbox{iff}(a_{i})=a_{i}\Leftrightarrow G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ ( $s\geq 0$ ). We have $I\models\mbox{iff}(a_{i})$ . So if ${\mathbf{u}}_{I}(i)=1$ , $I\models a_{i}$ , and hence $I\models G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ , giving $d_{i}\geq 1$ because $d_{i}$ is the number of rule bodies in $\{G_{j_{1}},\ldots,G_{j_{s}}\}$ that are true in $I$ . So $\mbox{min}_{1}(d_{i})=1$ holds. Otherwise if ${\mathbf{u}}_{I}(i)=0$ , we have $I\not\models a_{i}$ and $I\not\models G_{j_{1}}\vee\cdots\vee G_{j_{s}}$ . Consequently none of the rule bodies are true in $I$ and we have $d_{i}=\mbox{min}_{1}(d_{i})=0$ . Putting the two together, we have ${\mathbf{u}}_{I}(i)=d_{i}$ . Since $i$ is arbitrary, we conclude ${\mathbf{u}}_{I}={\rm min}_{1}({\mathbf{d}}_{I})$ , or equivalently $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})\|_{2}=0$ . The converse is similarly proved. Q.E.D.

Proposition 1 says that whether $I$ is a supported model of the program $P$ or not is determined by computing ${\mathbf{u}}_{I}-{\mbox{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 $P_{0}=(Q_{0},D_{0})$ in (5) having three rules $\{r_{1},r_{2},r_{3}\}$ , take a model $I_{0}=\{p,q\}$ over the ordered atom set ${\cal 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{u}}_{I_{0}}=[\,1\;\;1\;\;0\;]^{T}$ , ${\mathbf{u}}_{I_{0}}^{\delta}=[\,1\;\;1\;\;0\;\;0\;\;0\;\;1\,]^{T}$ , $1-{\mathbf{u}}_{I_{0}}^{\delta}=[\,0\;\;0\;\;1\;\;1\;\;1\;\;0\,]^{T}$ and finally $Q_{0}(1-{\mathbf{u}}_{I_{0}}^{\delta})=[\,0\;\;1\;\;0\,]^{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 $\mbox{min}_{1}(Q_{0}(1-{\mathbf{u}}_{I_{0}}^{\delta}))=[\,0\;\;1\;\;0\,]^{T}$ indicates that only the second rule body is false and the other two bodies are true in $I_{0}$ . So its bit inversion $1-\mbox{min}_{1}(Q_{0}(1-{\mathbf{u}}_{I_{0}}^{\delta}))=[\,1\;\;0\;\;1\,]^{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}}=D_{0}(1-\mbox{min}_{1}(Q_{0}(1-{\mathbf{u}}_{I_{0}}^{\delta})))=[\,1\;\;1\;\;0\,]^{T}$ .

${\mathbf{d}}_{I_{0}}=[\,1\;\;1\;\;0\,]^{T}$ denotes for each atom $a\in{\cal 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 ${\cal A}_{0}$ has one rule ( $p\leftarrow q\wedge\neg r$ ) whose body ( $q\wedge\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 $\mbox{min}_{1}({\mathbf{d}}_{I_{0}})=[\,1\;\;1\;\;0\,]^{T}$ denotes the truth values of the right hand sides of the completed rules $\{\mbox{iff}(p),\mbox{iff}(q),\mbox{iff}(r)\}$ evaluated by $I_{0}$ . Since ${\mathbf{u}}_{I_{0}}=\mbox{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},\ldots,r_{m}\}$ in a set ${\cal A}=\{a_{1},\ldots,a_{n}\}$ of $n$ ordered atoms as before. We first show the evaluation of the reduct of the matricized program $P=(D,Q)$ by a vectorized model ${\mathbf{u}}_{I}$ . Write $Q\in\{0,1\}^{m\times 2n}$ as $Q=[Q^{(1)}\;Q^{(2)}]$ where $Q^{(1)}\in\{0,1\}^{m\times n}$ (resp. $Q^{(2)}\in\{0,1\}^{m\times n}$ ) is the left half (resp. the right half) of $Q$ representing the positive literals (resp. negative literals) of each rule body in $Q$ . Compute $M^{(2)}=1-\mbox{min}_{1}(Q^{(2)}{\mathbf{u}}_{I})$ . It is an $m\times 1$ matrix (treated as a column vector here) such that $M^{(2)}(j)=0$ if the body of $r_{j}$ contains a negative literal false in $I$ and $M^{(2)}(j)=1$ otherwise ( $1\leq j\leq m$ ). Let $r_{j}^{+}$ be a rule $r_{j}$ with negative literals in the body deleted. We see that $P^{I}=\{r_{j}^{+}\mid M^{(2)}(j)=1,1\leq j\leq m\}$ and $P^{I}$ is syntactically represented by $(D^{I},Q^{(1)})$ where $D^{I}=D$ with columns $D(:,j)$ replaced by the zero column vector if $M^{(2)}(j)=0$ ( $1\leq j\leq m$ ). $D^{I}(i,:)$ denotes a rule set $\{r_{j}^{+}\mid D^{I}(i,j)=1,1\leq j\leq m\}$ in $P^{I}$ for $a_{i}\in{\cal A}$ . We call $P^{I}=(D^{I},Q^{(1)})$ the matricized reduct of $P$ by $I$ .

The matricized reduct $P^{I}=(D^{I},Q^{(1)})$ is evaluated in vector spaces as follows. Compute $M^{(1)}=M^{(2)}\odot(1-\mbox{min}_{1}(Q^{(1)}(1-{\mathbf{u}}_{I})))$ ⁷⁷7 $\odot$ is component-wise product. . $M^{(1)}$ denotes the truth values of rule bodies in $P^{I}$ evaluated by $I$ . Thus $M^{(1)}(j)=1$ ( $1\leq j\leq m$ ) if $r_{j}^{+}$ is contained in $P^{I}$ and its body is true in $I$ . Otherwise $M^{(1)}(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}^{+}=DM^{(1)}$ . ${\mathbf{d}}_{I}^{+}(i)$ ( $1\leq i\leq n$ ) is the number of rules in $P^{I}$ for the $i$ -th atom $a_{i}$ in ${\cal A}$ whose bodies are true in $I$ .

Proposition 2

Let $P=(D,Q)$ be a matricized program $P$ in a set ${\cal A}=\{a_{1},\ldots,a_{n}\}$ of $n$ ordered atoms and $I$ a model over ${\cal A}$ . Write $Q=[Q^{(1)}\;Q^{(2)}]$ as above. Let ${\mathbf{u}}_{I}$ be the vectorized model $I$ . Compute $M^{(2)}=1-{\rm min}_{1}(Q^{(2)}{\mathbf{u}}_{I})$ , $M^{(1)}=M^{(2)}\odot(1-{\rm min}_{1}(Q^{(1)}(1-{\mathbf{u}}_{I})))$ and ${\mathbf{d}}_{I}^{+}=DM^{(1)}$ . Also compute ${\mathbf{d}}_{I}=D(1-{\rm min}_{1}(Q(1-{\mathbf{u}}_{I}^{\delta})))$ . Then, $I\models{\rm comp}(P)$ , $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})\|_{2}=0$ , $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I}^{+})\|_{2}=0$ and $I\models{\rm comp}(P^{I})$ are all equivalent.

(Proof) We prove ${\mathbf{d}}_{I}={\mathbf{d}}_{I}^{+}$ first. Recall that a rule $r_{j}^{+}$ in $P^{I}$ is created by removing negative literals true in $I$ from the body of $r_{j}$ in $P$ . So for any $a_{i}\in{\cal A}$ , it is immediate that $a_{i}$ has a rule $r_{j}\in P$ whose body is true in $I$ iff $a_{i}$ has the rule $r_{j}^{+}\in P^{I}$ whose body is true in $I$ . Thus ${\mathbf{d}}_{I}(i)={\mathbf{d}}_{I}^{+}(i)$ for every $i$ ( $1\leq i\leq n$ ), and hence ${\mathbf{d}}_{I}={\mathbf{d}}_{I}^{+}$ . Consequently, we have $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})\|_{2}=0$ iff $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I}^{+})\|_{2}=0$ . Also $I\models\mbox{comp}(P^{I})$ iff $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I}^{+})\|_{2}=0$ is proved similarly to Proposition 1 (details omitted). The rest follows from Proposition 1. Q.E.D.

From the viewpoint of end-to-end ASP, Proposition 2 means that we can obtain a supported model $I$ as a binary solution ${\mathbf{u}}_{I}$ of the equation ${\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})=0$ derived from $P$ or ${\mathbf{u}}_{I}-{\rm 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{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})$ is piecewise linear w.r.t. ${\mathbf{u}}_{I}$ whereas the latter ${\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I}^{+})$ is piecewise quadratic w.r.t. ${\mathbf{u}}_{I}$ .

Now look at $P_{0}=\{r_{1},r_{2},r_{3}\}$ in (5) and a model $I_{0}=\{p,q\}$ again. $\begin{array}[]{lcl}P_{0}^{I_{0}}&=&\begin{array}[]{ll}\left\{\;\begin{array}[]{lllll}p&\leftarrow&q\\ q&\leftarrow&\\ \end{array}\right.\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 $P_{0}=(D_{0},Q_{0})$ be the matricized $P_{0}$ . We first decompose $Q_{0}$ in (8) as $Q_{0}=[Q_{0}^{(1)}\;Q_{0}^{(2)}]$ where $Q_{0}^{(1)}$ is the positive part and $Q_{0}^{(2)}$ the negative part of $Q_{0}$ . They are

\displaystyle\hskip 20.00003pt\begin{array}[]{lcl}Q_{0}^{(1)}&=&\begin{bmatrix}\;\;0&\,1&\,0\;\;\\ \;\;0&\,0&\,0\;\;\\ \;\;0&\,0&\,0\;\;\\ \end{bmatrix}\end{array}\hskip 10.00002pt\mbox{and}\hskip 10.00002pt\begin{array}[]{lcl}Q_{0}^{(2)}&=&\begin{bmatrix}\;\;0&\,0&\,1\;\;\\ \;\;0&\,1&\,0\;\;\\ \;\;0&\,0&\,0\;\;\\ \end{bmatrix}\end{array}.

Let ${\mathbf{u}}_{I_{0}}=[\,1\;\;1\;\;0\,]^{T}$ be the vectorized $I_{0}$ . We first compute $M_{0}^{(2)}=1-\mbox{min}_{1}(Q_{0}^{(2)}{\mathbf{u}}_{I_{0}})$ to determine rules to be removed. Since $M_{0}^{(2)}=[\,1\;\;0\;\;1\,]^{T}$ , the second rule $r_{2}$ , indicated by $M_{0}^{(2)}(2)=0$ , is removed from $P_{0}$ , giving $P_{0}^{I_{0}}=\{r_{1}^{+},r_{3}^{+}\}$ . Using $M_{0}^{(2)}$ and $D_{0}$ shown in (11), we then compute $M_{0}^{(1)}=M_{0}^{(2)}\odot(1-\mbox{min}_{1}(Q_{0}^{(1)}(1-{\mathbf{u}}_{I_{0}})))=[\,1\;\;0\;\;1\,]^{T}$ and ${\mathbf{d}}_{I_{0}}^{+}=D_{0}M_{0}^{(1)}=[\,1\;\;1\;\;0\,]^{T}$ . ${\mathbf{d}}_{I_{0}}^{+}$ denotes the number rule bodies in $P_{0}^{I_{0}}$ true in $I_{0}$ for each atom. Thus, since ${\mathbf{u}}_{I_{0}}=\mbox{min}_{1}({\mathbf{d}}_{I_{0}}^{+})(=[\,1\;\;1\;\;0\,]^{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 [2]. They encode normal logic programs differently from ours based on Sakama’s encoding [22] 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 [27] which encodes a program in terms of a single matrix and evaluates conjunctions based on the number of true literals. They compute supported models by minimizing a non-negative function, not solving equations like [2]. 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{u}}_{I}$ representing a supported model $I$ of a matricized program $P=(D,Q)$ that satisfies $\|{\mathbf{u}}_{I}-{\rm min}_{1}({\mathbf{d}}_{I})\|_{2}=0$ where ${\mathbf{d}}_{I}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}_{I}^{\delta})))$ . For this task, we relax ${\mathbf{u}}_{I}\in\{1,0\}^{n}$ to ${\mathbf{u}}\in\mathbb{R}^{n}$ and introduce a non-negative cost function $J_{SU}$ :

	$\displaystyle J_{SU}$	$\displaystyle=$	$\displaystyle 0.5\cdot\left(\\|{\mathbf{u}}-{\rm min}_{1}({\mathbf{d}})\\|_{2}^{2}+\ell_{2}\cdot\\|{\mathbf{u}}\odot(1-{\mathbf{u}})\\|_{2}^{2}\right)$
			$\displaystyle\;\mbox{where}\;\;\ \ell_{2}>0\;\mbox{and}\;{\mathbf{d}}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}^{\delta}))).$

Proposition 3

Let $J_{SU}$ be defined from a program $P=(D,Q)$ as above.

$J_{SU}=0$ if-and-only-if ${\mathbf{u}}$ is a binary vector representing a supported model of $P$ .

(Proof) Apparently if $J_{SU}=0$ , we have $\|{\mathbf{u}}-{\rm min}_{1}({\mathbf{d}})\|_{2}^{2}=0$ and $\|{\mathbf{u}}\odot(1-{\mathbf{u}})\|_{2}^{2}=0$ . The second equation means ${\mathbf{u}}$ is binary ( $x(1-x)=0\Leftrightarrow x\in\{1,0\}$ ), and the first equation means this binary ${\mathbf{u}}$ is a vector representing a supported model of $P$ by Proposition 1. The converse is obvious. Q.E.D.

$J_{SU}$ is piecewise differentiable and we can obtain a supported model of $P$ as a root ${\mathbf{u}}$ of $J_{SU}$ by minimizing $J_{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},\ldots,a_{n}\}$ and ${\mathbf{u}}=[u_{1},\ldots,u_{n}]^{T}$ represents their continuous truth values where ${\mathbf{u}}(p)=u_{p}\in\mathbb{R}$ is the continuous truth value for atom $a_{p}$ ( $1\leq p\leq n$ ). For the convenience of derivation, we introduce the dot product $({A}\bullet{B})=\sum_{i,j}{A}(i,j){B}(i,j)$ of matrices $A$ and $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 $({A}\bullet({B}\odot{C}))=(({B}\odot{A})\bullet{C})$ and $({A}\bullet({B}{C}))=(({B}^{T}{A})\bullet{C})=(({A}{C}^{T})\bullet{B})$ hold.

Let $P=(D,Q)$ be the matricized program and write $Q=[Q^{(1)}\,Q^{(2)}]$ . Introduce $N$ , $M$ , ${\mathbf{d}}$ , $E$ , $F$ and compute $J_{SU}$ by

\displaystyle\begin{array}[]{lclll}N&=&Q(1-{\mathbf{u}}^{\delta})=Q^{(1)}(1-{\mathbf{u}})+Q^{(2)}{\mathbf{u}}\\ &&\mbox{: (continuous) counts of false literals in the rule bodies}\\ M&=&1-{\rm min}_{1}(N)\\ &&\mbox{: (continuous) truth values of the rule bodies}\\ {\mathbf{d}}&=&DM\\ &&\mbox{: (continuous) counts of true disjuncts for each atom}\\ E&=&{\rm min}_{1}({\mathbf{d}})-{\mathbf{u}}\\ &&\mbox{: error between the estimated truth values of atoms and ${\mathbf{u}}$}\\ F&=&{\mathbf{u}}\odot(1-{\mathbf{u}})\\ J_{sq}&=&(E\bullet E)\\ J_{nrm}&=&(F\bullet F)\\ J_{SU}&=&J_{sq}+\ell_{2}\cdot J_{nrm}.\end{array}

(27)

We then compute the Jacobian $J_{a_{SU}}$ of $J_{SU}$ . We first compute $\displaystyle{\frac{\partial J_{sq}}{\partial{u}_{p}}}$ where ${u}_{p}={\mathbf{u}}(p)$ ( $1\leq p\leq n$ ).

$\displaystyle\frac{\partial M}{\partial u_{p}}$	$\displaystyle=$	$\displaystyle-[N\leq 1]\odot\bigl{(}(Q^{(2)}-Q^{(1)}){\mathbf{I}}_{p}\bigr{)}=[N\leq 1]\odot\bigl{(}(Q^{(1)}-Q^{(2)}){\mathbf{I}}_{p}\bigr{)}$
$\displaystyle\frac{\partial J_{sq}}{\partial u_{p}}$	$\displaystyle=$	$\displaystyle\Bigl{(}E\,\bullet\;[DM\leq 1]\odot\Bigl{(}D\left(\frac{\partial M}{\partial u_{p}}\right)\Bigr{)}-{\mathbf{I}}_{p}\Bigr{)}$
	$\displaystyle=$	$\displaystyle\Bigl{(}E\,\bullet\;[DM\leq 1]\odot(D([N\leq 1]\odot((Q^{(1)}-Q^{(2)}){\mathbf{I}}_{p})))-{\mathbf{I}}_{p}\Bigr{)}$
	$\displaystyle=$	$\displaystyle\Bigl{(}D^{T}([DM\leq 1]\odot E)\;\bullet\;[N\leq 1]\odot(((Q^{(1)}-Q^{(2)}){\mathbf{I}}_{p})\Bigr{)}-(E\;\bullet\;{\mathbf{I}}_{p})$
	$\displaystyle=$	$\displaystyle\Bigl{(}(Q^{(1)}-Q^{(2)})^{T}([N\leq 1]\odot(D^{T}([DM\leq 1]\odot E)))-E\;\bullet\;{\mathbf{I}}_{p}\Bigr{)}$

Since $p$ is arbitrary, we have ${\displaystyle\frac{\partial J_{sq}}{\partial{\mathbf{u}}}=(Q^{(1)}-Q^{(2)})^{T}([N\leq 1]\odot(D^{T}([DM\leq 1]\odot E)))-E}$ .

Next we compute ${\displaystyle\frac{\partial J_{nrm}}{\partial u_{p}}}$ :

$\displaystyle\frac{\partial F}{\partial u_{p}}$	$\displaystyle=$	$\displaystyle\left(\frac{\partial{\mathbf{u}}}{\partial u_{p}}\right)\odot(1-{\mathbf{u}})+{\mathbf{u}}\odot\left(\frac{\partial(1-{\mathbf{u}})}{\partial u_{p}}\right)$
	$\displaystyle=$	$\displaystyle({\mathbf{I}}_{p}\odot(1-{\mathbf{u}}))-({\mathbf{u}}\odot{\mathbf{I}}_{p})=(1-2{\mathbf{u}})\odot{\mathbf{I}}_{p}$
$\displaystyle\frac{\partial J_{nrm}}{\partial u_{p}}$	$\displaystyle=$	$\displaystyle\bigl{(}F\;\bullet\;\left(\frac{\partial F}{\partial u_{p}}\right)\bigr{)}$
	$\displaystyle=$	$\displaystyle(F\;\bullet\;(1-2{\mathbf{u}})\odot{\mathbf{I}}_{p})=((1-2{\mathbf{u}})\odot F\;\bullet\;{\mathbf{I}}_{p})$

Again since $p$ is arbitrary, we have ${\displaystyle\frac{\partial J_{nrm}}{\partial{\mathbf{u}}}=(1-2{\mathbf{u}})\odot F}$ and reach

$\displaystyle J_{a_{SU}}$	$\displaystyle=$	$\displaystyle\left(\frac{\partial J_{sq}}{\partial{\mathbf{u}}}\right)+\ell_{2}\cdot\left(\frac{\partial J_{nrm}}{\partial{\mathbf{u}}}\right)$
	$\displaystyle=$	$\displaystyle(Q^{(1)}-Q^{(2)})^{T}([N\leq 1]\odot(D^{T}([{\mathbf{d}}\leq 1]\odot E)))-E+\ell_{2}\cdot((1-2{\mathbf{u}})\odot F)$
		$\displaystyle\;\mbox{where}\;\;N=Q(1-{\mathbf{u}}^{\delta}),\,{\mathbf{d}}=D(1-{\rm min}_{1}(N)),\,E={\rm min}_{1}({\mathbf{d}})-{\mathbf{u}},\,$
		$\displaystyle\;\;\mbox{and}\;\;F={\mathbf{u}}\odot(1-{\mathbf{u}}).$

3.4 Adding constraints

A rule which has no head like $\leftarrow a\wedge\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 just like rule bodies in a program using a binary matrix $Q_{c}=[Q_{c}^{(1)}Q_{c}^{(2)}]$ . We call $Q_{c}$ constraint matrix. We introduce $N_{c}$ , a non-negative function $J_{c}$ of ${\mathbf{u}}$ and $J_{c}$ ’s Jacobian $J_{a_{c}}$ as follows.

$\displaystyle N_{c}$	$\displaystyle=$	$\displaystyle Q_{c}(1-{\mathbf{u}}^{\delta})=Q_{c}^{(1)}(1-{\mathbf{u}})+Q_{c}^{(2)}{\mathbf{u}}$
$\displaystyle J_{c}$	$\displaystyle=$	$\displaystyle({\mathbf{1}}\bullet(1-{\rm min}_{1}(N_{c})))\hskip 5.0pt\mbox{where ${\mathbf{1}}$ is an all-ones vector}$	(29)
$\displaystyle J_{a_{c}}$	$\displaystyle=$	$\displaystyle(Q_{c}^{(1)}-Q_{c}^{(2)})^{T}[N_{c}\leq 1]$	(30)

The meaning of $N_{c}$ and $J_{c}$ is clear when ${\mathbf{u}}$ is binary. Note that any binary ${\mathbf{u}}$ is considered as a model over a set ${\cal A}=\{a_{1},\ldots,a_{n}\}$ of $n$ ordered atoms in an obvious way. Suppose $k$ constraints are given to be satisfied. Then $Q_{c}$ is a $k\times 2n$ binary matrix and $N_{c}$ is a $k\times 1$ matrix. $N_{c}(i)$ ( $1\leq i\leq k$ ) is the number of literals falsified by ${\mathbf{u}}$ in a conjunction $G_{i}$ of the $i$ -th constraint $\leftarrow G_{i}$ . So $N_{c}(i)=0$ , or equivalently $1-{\rm min}_{1}(N_{c}(i))=1$ implies $G_{i}$ has no false literal i.e. ${\mathbf{u}}\models G_{i}$ , and vice versa. Hence $J_{c}=\sum_{i=1}^{k}(1-{\rm min}_{1}(N_{c}(i)))=({\mathbf{1}}\bullet(1-{\rm min}_{1}(N_{c})))$ equals the number of violated constraints. Consequently when ${\mathbf{u}}$ is binary. we can say that $J_{c}=0$ iff all constraints are satisfied by ${\mathbf{u}}$ .

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

3.5 An algorithm for computing supported models with constraints

Here we present a minimization algorithm for computing supported models of the matricized program $P=(D,Q)$ which satisfy constraints represented by a constraint matrix $Q_{c}$ . We first combine $J_{SU}$ and $J_{c}$ into $J_{SU+c}$ .

$\displaystyle J_{SU+c}$	$\displaystyle=$	$\displaystyle J_{SU}+\ell_{3}\cdot J_{c}$
	$\displaystyle=$	$\displaystyle 0.5\cdot\left(\\|{\mathbf{u}}-{\rm min}_{1}({\mathbf{d}})\\|_{2}^{2}+\ell_{2}\cdot\\|{\mathbf{u}}\odot(1-{\mathbf{u}})\\|_{2}^{2}\right)$
		$\displaystyle+\ell_{3}\cdot({\mathbf{1}}\bullet(1-{\rm min}_{1}(Q_{c}(1-{\mathbf{u}}^{\delta}))))\hskip 10.00002pt\ell_{2}>0,\ell_{3}>0$
		$\displaystyle\;\mbox{where}\;\;{\mathbf{d}}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}^{\delta})))$
$\displaystyle J_{a_{SU+c}}$	$\displaystyle=$	$\displaystyle J_{a_{SU}}+\ell_{3}\cdot J_{a_{c}}$	(32)

The next proposition is immediate from Proposition 3.

Proposition 4

$J_{SU+c}=0$ if-and-only-if ${\mathbf{u}}$ represents a supported model of $P$ satisfying a constraint matrix $Q_{c}$ .

We compute $J_{SU}$ in $J_{SU+c}$ by (27) and $J_{c}$ by (29), and their Jacobians $J_{a_{SU}}$ and $J_{a_{c}}$ by (3.3) and by (30) respectively. We minimize the non-negative $J_{SU+c}$ to zero by Newton’s method using Algorithm 1. It finds a solution ${\mathbf{u}}_{*}$ of $J_{SU+c}=0$ which represents a supported model of $P$ satisfying constraint matrix $Q_{c}$ . The updating formula is derived from the first order Taylor expansion of $J_{SU+c}$ and by solving $J_{SU+c}+(J_{a_{SU*c}}\bullet({\mathbf{u}}_{\rm new}-{\mathbf{u}}))=0$ w.r.t. ${\mathbf{u}}_{\rm new}$ . We use it with a learning rate $\alpha>0$ as follows.

\displaystyle{\mathbf{u}}_{\rm new}

\displaystyle=

\displaystyle{\mathbf{u}}-\alpha\left(\frac{J_{SU+c}}{(J_{a_{SU+c}}\;\bullet\;J_{a_{SU+c}})}\right)J_{a_{SU+c}}

(33)

1 Input: matricized program

P=(D,Q)

, constraint matrix

Q_{c}

max\_itr\in\mathbb{Z}

max\_try\in\mathbb{Z}

2 Output: binary vector

{\mathbf{u}}_{*}

representing a supported model of

P

satisfying constraints represented by

Q_{c}

{\mathbf{u}}\leftarrow

random initialization

4 for $i\leftarrow 1$ to max_try do

5 for $j\leftarrow 1$ to max_itr do

6 optimally threshold

{\mathbf{u}}

to a binary vector

{\mathbf{u}}_{*}

so that

error\leftarrow\|{\mathbf{u}}_{*}-{\rm min}_{1}({\mathbf{d}}_{*})\|_{2}^{2}+({\mathbf{1}}\bullet(1-{\rm min}_{1}(Q_{c}(1-{\mathbf{u}}_{*}^{\delta}))))

8 is minimum where

{\mathbf{d}}_{*}=D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}_{*}^{\delta}))

9 if error = 0 then

break

10Update

{\mathbf{u}}

by (33)

11if error = 0 then

break

12perturbate

{\mathbf{u}}

to escape from a local minimum

13return ${\mathbf{u}}_{*}$

Algorithm 1 minimizing

J_{SU+c}

to zero

Algorithm 1 is a double loop algorithm where the inner $j$ -loop updates ${\mathbf{u}}\in\mathbb{R}^{n}$ repeatedly to minimize $J_{SU+c}$ while thresholding ${\mathbf{u}}$ into a binary solution candidate ${\mathbf{u}}_{*}\in\{1,0\}^{n}$ for $J_{SU+c}=0$ . The outer $i$ -loop is for retry when the inner fails to find a solution. The initialization at line 3 is carried out by sampling ${\mathbf{u}}(i)\sim{\cal N}(0,1)+0.5$ ( $1\leq i\leq n$ ) where ${\cal N}(0,1)$ is the standard normal distribution. Lines 6,7 and 8 collectively perform thresholding of ${\mathbf{u}}$ into a binary ${\mathbf{u}}_{*}$ . As the inner loop repeats, $J_{SU+c}$ becomes smaller and smaller and so do $J_{sq}$ and $J_{nrm}$ in $J_{SU}$ . $J_{sq}$ being small means ${\mathbf{u}}$ is close to a supported model of $P$ while $J_{nrm}$ being small means each element of ${\mathbf{u}}$ is close to $\{1,0\}$ . So binarization ${\mathbf{u}}_{*}=[{\mathbf{u}}\geq\theta]$ with an appropriate threshold $\theta$ ⁸⁸8 Currently given ${\mathbf{u}}$ , we divide the interval $[{\rm min}({\mathbf{u}}),{\rm max}({\mathbf{u}})]$ into 20 equally distributed notches and use each notch as a threshold value $\theta$ . has a good chance of yielding a binary ${\mathbf{u}}_{*}$ representing a supported model of $P$ satisfying constraints represented by $Q_{c}$ . It may happen that the inner loop fails to find a solution. In such case, we retry another $j$ -loop with perturbated ${\mathbf{u}}$ at line 12. There ${\mathbf{u}}$ is perturbated by ${\mathbf{u}}\leftarrow 0.5({\mathbf{u}}+\Delta+0.5)$ where $\Delta\sim{\cal N}(0,1)$ before the next $j$ -loop.

3.6 Connection to neural network computation

At this point, it is quite interesting to see the connection of our approach to neural network computation. In (27), we compute $M$ and ${\mathbf{d}}=DM$ . 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{u}}$ . Consider the computation of $M=(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}^{\delta})))$ . We rewrite this using $1-{\rm min}(x,1)={\rm ReLU}(1-x)$ to

$\displaystyle M$	$\displaystyle=$	$\displaystyle 1-\mbox{min}_{1}(Q^{(1)}(1-{\mathbf{u}})+Q^{(2)}{\mathbf{u}})$
	$\displaystyle=$	$\displaystyle{\rm ReLU}(W{\mathbf{u}}+{\mathbf{b}})$
		$\displaystyle\mbox{where}\;\;Q=[Q^{(1)}\,Q^{(2)}],\,W=Q^{(1)}-Q^{(2)},\,{\mathbf{b}}=1-Q^{(1)}{\mathbf{1}}$

So $M$ is the output of a ReLU network having a weight matrix $W=Q^{(1)}-Q^{(2)}$ and a bias vector ${\mathbf{b}}=1-Q^{(1)}{\mathbf{1}}$ . Then ${\rm min}_{1}({\mathbf{d}})={\rm min}_{1}(DM)={\rm min}_{1}(D\cdot{\rm ReLU}(W{\mathbf{u}}+{\mathbf{b}}))$ is the output of a ReLU network with a single hidden layer and a linear output layer represented by $D$ having ${\rm min}_{1}(\cdot)$ as activation function. Also when we compute a stable model ${\mathbf{u}}$ , we minimize $J_{SU+c}$ (27) which contains an MSE error term $J_{sq}=\|{\rm min}_{1}({\mathbf{d}})-{\mathbf{u}}\|^{2}$ using $J_{a_{SU+c}}$ (32). This is precisely back propagation from learning data ${\mathbf{u}}$ .

Thus we may say that our approach is an integration of ASP semantics and neural computation and provides a neuro-symbolic [23] way of ASP computation. Nonetheless, there is a big difference. In standard neural network architecture, a weight matrix $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.

4 Computing stable models in vector spaces

4.1 Loop formulas and stable models

Let $P=(D,Q)$ be a matricized program in a set of atoms ${\cal A}=\{a_{1},\ldots,a_{n}\}$ having $m$ rules $\{a_{i_{1}}\leftarrow G_{1},\ldots,a_{i_{m}}\leftarrow G_{m}\}$ where $D\in\{1,0\}^{n\times m}$ and $Q\in\{1,0\}^{m\times 2n}$ . We assume atoms and rules are ordered as indicated.

Computing a supported model of $P$ is equivalent to computing any binary fixedpoint ${\mathbf{u}}\in\{1,0\}^{n}$ such that ${\mathbf{u}}={\rm min}_{1}(D(1-\mbox{min}_{1}(Q(1-{\mathbf{u}}^{\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 [14] 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 $L=\{h_{1},\ldots,h_{k}\}\subseteq{\cal A}$ be a loop in $P$ . Recall that $L$ is a set of atoms which are strongly connected in the positive dependency graph of $P$ ⁹⁹9 In the case of a singleton loop $L=\{h\}$ , we specifically require, following [14], that $h$ has a self-loop, i.e. there must be a rule of the form $h\leftarrow h\wedge H$ in $P$ . . A support rule for $h$ is a rule $h\leftarrow H$ such that $H^{+}\cap L=\emptyset$ . $H$ is called a support body for $L$ . Introduce a loop formula for $L$ ¹⁰¹⁰10 In the original form, the antecedent of $LF(L)$ is a disjunction $(h_{1}\vee\cdots\vee h_{p})$ [14]. Later it is shown that Lin-Zhao’s theorem also holds for the AND type $LF(L)$ [6]. We choose this AND type $LF(L)$ as it is easier to satisfy. by

	$\displaystyle LF(L)$	$\displaystyle=$	$\displaystyle(h_{1}\wedge\cdots\wedge h_{p})\rightarrow(H_{1}\vee\cdots\vee H_{q})$		(34)
			$\displaystyle\mbox{where $\{H_{1},\ldots,H_{q}\}$ are support bodies for $L$}.$		(34)

Then define loop formulas associated with $P$ as $LF(P)=\{LF(L)\mid\mbox{$L$ is a loop in $P$}\}$ . Logically $LF(P)$ is treated as the conjunction of its elements and we sometimes call it the loop formula associated with $P$ . Now we evaluate $LF(P)$ by a real vector ${\mathbf{u}}\in\mathbb{R}^{n}$ . Introduce an external support matrix $E_{s}\in\{1,0\}^{n\times m}$ by $E_{s}(i,j)=1$ if there is a support rule $a_{i}\leftarrow G_{j}$ for $a_{i}\in{\cal A}$ , else $E_{s}(i,j)=0$ ( $1\leq i\leq n,1\leq j\leq m$ ). Suppose there are $t$ loops $\{L_{1},\ldots,L_{t}\}$ in $P$ . Introduce a loop matrix $L_{oop}\in\{1,0\}^{t\times m}$ such that $L_{oop}(s,j)=1$ if the $s$ -th loop $L_{s}$ has $G_{j}$ as a support body for $L_{s}$ , else $L_{oop}(s,j)=0$ ( $1\leq s\leq t$ ). Evaluate $J_{LF}$ by ${\mathbf{u}}$ as follows.

$\displaystyle M$	$\displaystyle=\;$	$\displaystyle 1-{\rm min}_{1}(Q(1-{\mathbf{u}}^{\delta}))\hskip 20.00003pt\mbox{: (continuous) truth values by ${\mathbf{u}}$ of the rule bodies in $P$}$
$\displaystyle L_{s}$	$\displaystyle=\;$	$\displaystyle L_{oop}(s,:)\hskip 70.0001pt\mbox{: represents the $s$-th loop in $\{L_{1},\ldots,L_{t}\}$}$
$\displaystyle A_{s}$	$\displaystyle=\;$	$\displaystyle L_{s}(1-{\mathbf{u}})+L_{s}E_{s}M\hskip 30.00005pt\mbox{: (continuous) counts of true disjuncts by ${\mathbf{u}}$ of $LF(L_{s})$ }$
$\displaystyle J_{LF}$	$\displaystyle=\;$	$\displaystyle\sum_{s=1}^{t}\left(1-{\rm min}_{1}(A_{s})\right)$	(35)

$J_{LF}$ is a non-negative piecewise linear function of ${\mathbf{u}}$ .

Proposition 5

Let $J_{LF}$ be defined as above. When ${\mathbf{u}}$ is binary, it holds that

$J_{LF}=0$ if-and-only-if ${\mathbf{u}}\models LF(P)$ .

(Proof) Suppose $J_{LF}=0$ and ${\mathbf{u}}$ is binary. A summand $(1-{\rm min}_{1}(A_{s}))$ in $J_{LF}$ (35) corresponds to the $s$ -th loop $L_{s}=\{h_{1},\ldots,h_{p}\}$ and is non-negative. Consider $LF(L_{s})=(h_{1}\wedge\cdots\wedge h_{p})\rightarrow(H_{1}\vee\cdots\vee H_{q})$ as a disjunction $\neg h_{1}\vee\cdots\vee\neg h_{p}\vee H_{1}\vee\cdots\vee H_{q}$ . Then $J_{LF}=0$ implies $(1-{\rm min}_{1}(A_{s}))=0$ , or equivalently $A_{s}\geq 1$ . Consequently, as ${\mathbf{u}}$ is binary, we have $L_{s}(1-{\mathbf{u}})\geq 1$ or $L_{s}E_{s}M\geq 1$ . The former means ${\mathbf{u}}\models\neg h_{1}\vee\cdots\vee\neg h_{p}$ . The latter, $L_{s}E_{s}M\geq 1$ , means ${\mathbf{u}}\models H_{1}\vee\cdots\vee H_{q}$ . This is because $(E_{s}M)(i)$ is the number of support rules for $a_{i}\in{\cal A}$ whose bodies are true in ${\mathbf{u}}$ ( $1\leq i\leq n$ ), and hence $L_{s}E_{s}M\geq 1$ means some support body $H_{r}$ ( $1\leq r\leq q$ ) for $L_{s}$ is true in ${\mathbf{u}}$ . So in either case ${\mathbf{u}}\models LF(L_{s})$ . Since $s$ is arbitrary, we have ${\mathbf{u}}\models LF(P)$ . The converse is straightforward and omitted. Q.E.D.

The Jacobian $J_{a_{LF}}$ of $J_{LF}$ is computed as follows.

$\displaystyle N$	$\displaystyle=$	$\displaystyle Q(1-{\mathbf{u}}^{\delta})$
$\displaystyle N_{s}$	$\displaystyle=$	$\displaystyle L_{s}(1-{\mathbf{u}})$
$\displaystyle M_{s}$	$\displaystyle=$	$\displaystyle{\rm min}_{1}(N_{s})$
$\displaystyle J_{a_{LF}}$	$\displaystyle=$	$\displaystyle\frac{\partial J_{LF}}{\partial{\mathbf{u}}}=\sum_{s=1}^{t}-\left(\frac{\partial{\rm min}_{1}(A_{s})}{\partial{\mathbf{u}}}\right)$	(36)
	$\displaystyle=$	$\displaystyle-\sum_{s=1}^{t}[A_{s}\leq 1]\left(\left(\frac{\partial M_{s}}{\partial{\mathbf{u}}}\right)+L_{s}E_{s}\left(\frac{\partial M}{\partial{\mathbf{u}}}\right)^{T}\right)$
	$\displaystyle=$	$\displaystyle\sum_{s=1}^{t}[A_{s}\leq 1]\left([N_{s}\leq 1]L_{s}^{T}+(((L_{s}E_{s})\odot[N\leq 1]^{T})(Q^{(2)}-Q^{(1)})^{T}\right)$

Here $Q=[Q^{(1)}\,Q^{(2)}]$ and $L_{s}$ , $A_{s}$ and $M$ are computed by (35).

Now introduce a new cost function $J_{SU+c+LF}$ by (37) that incorporates $J_{LF}$ and compute its Jacobian $J_{a_{SU+c+LF}}$ by (32).

	$\displaystyle J_{SU+c+LF}$	$\displaystyle=$	$\displaystyle J_{SU+c}+\ell_{4}\cdot J_{LF}\;\;\mbox{where $\ell_{4}>0$}$		(37)
	$\displaystyle J_{a_{SU+c+LF}}$	$\displaystyle=$	$\displaystyle J_{a_{SU+c}}+\ell_{4}\cdot J_{a_{LF}}$		(38)

By combining Proposition 4, 5 and Lin-Zhao’s theorem [14], the following is obvious.

Proposition 6

${\mathbf{u}}$ is a stable model of $P$ satisfying constraints represented by $Q_{c}$ if-and-only-if ${\mathbf{u}}$ is a root of $J_{SU+c+LF}$ .

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

When a program $P$ is tight [5], for example when rules have no positive literal in their bodies, $P$ has no loop and hence $LF$ is empty. In such case, we directly minimize $J_{SU+c}$ instead of using $J_{SU+c+LF}$ with the empty $LF$ .

4.2 LF heuristics

Minimizing $J_{SU+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 [12] and they make the computation of $J_{LF}$ (35) extremely difficult or practically impossible. To mitigate this seemingly insurmountable difficulty, we propose two heuristics which use a subset of loop formulas.

$LF_{max}$ :

The first heuristics is $LF_{max}$ heuristics. We consider only a set $LF_{max}$ of loop formulas associated with SCCs in the positive dependency graph $\mbox{pdg}(P)=(V,E)$ of a program $P$ . In the case of a singleton SCC {a}, ’a’ must have a self-loop in $\mbox{pdg}(P)$ . We compute SCCs in $O(|E|+|V|)$ time by Tarjan’s algorithm [28].
$LF_{min}$ :

In this heuristics, instead of SCCs (maximal strongly connected subgraphs), we choose minimal strongly connected subgraphs, i.e. cycle graphs. Denote by $LF_{min}$ the set of loop formulas associated with cycle graphs in $\mbox{pdg}(P)$ . We use an enumeration algorithm described in [15] to enumerate cycles and construct $LF_{min}$ due to its simplicity.

We remark that although $LF_{max}$ and $LF_{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 $J_{SU+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 ${\cal 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^{+})$ , i.e. 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 ${\cal F}_{P}={\cal A}\setminus LM(P^{+})$ from our consideration as they are known to be false in any stable model. We call ${\cal 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^{+}$ , i.e. the total number of occurrences of atoms in $P^{+}$ [3]¹¹¹¹11 We implemented the linear time algorithm in [3] linear algebraically using vector and matrix and confirmed its linearity. . 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 $C$ , we can obtain downsized ones, $P^{\prime}$ and $C^{\prime}$ , as follows.

Step 1:

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

Step 2:

Define

$\displaystyle G^{\prime}$	$\displaystyle=$	conjunction $G$ with negative literals $\{\neg a\in G\mid a\in{\cal F}_{P}\}$ removed
$\displaystyle P^{\prime}$	$\displaystyle=$	$\displaystyle\{a\leftarrow G^{\prime}\mid a\leftarrow G\in P,a\not\in{\cal F}_{P},G^{+}\cap{\cal F}_{P}=\emptyset\}$
		$\displaystyle\;\mbox{where}\;\;G^{+}=\mbox{positive literals in $G$}$
$\displaystyle C^{\prime}$	$\displaystyle=$	$\displaystyle\{\leftarrow G^{\prime}\mid\leftarrow G\in C,G^{+}\cap{\cal F}_{P}=\emptyset\}$	(40)

Proposition 7

Let $P^{\prime}$ and $C^{\prime}$ be respectively the program (Step 2:) and constraints (40). 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 ${\cal F}_{P}$ is false in $I$ . Then

$I^{\prime}$ is a stable model of $P^{\prime}$ satisfying constraints $C^{\prime}$ if-and-only-if $I$ is a stable model of $P$ satisfying constraints $C$ .

(Proof) We prove first $I^{\prime}$ is a stable model of $P^{\prime}$ if-and-only-if $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\{\neg b\mid b\in{\cal F}_{P}\}$ and $G^{+}\cap{\cal F}_{P}=\emptyset$ . $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\{\neg b\mid b\in{\cal 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 an 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{\cal F}_{P}=\emptyset$ , i.e., $LM(P^{I})$ contains no stable false atom. So $b\not\in{\cal F}_{P}$ and $G^{+}\cap{\cal F}_{P}=\emptyset$ because every atom in the SDL 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}$ . We can also similarly prove that $I^{\prime}$ satisfies $C^{\prime}$ if-and-only-if $I$ satisfies $C$ . So we are done. Q.E.D.

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 paper, the program size is mostly small¹²¹²12 Matricized programs in this paper are all written in GNU Octave 6.4.0 and run on a PC with Intel(R) Core(TM) [email protected] CPU with 26GB memory. .

5.1 The 3-coloring problem

We first deal with the 3-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) in Fig. 1.

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}({\rm red}),a_{2}({\rm blue}),a_{3}({\rm 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.

\displaystyle P_{1}

\displaystyle=

\displaystyle\left\{\;\begin{array}[]{lllll}a_{1}\leftarrow\neg a_{2}\wedge\neg a_{3},&a_{2}\leftarrow\neg a_{3}\wedge\neg a_{1},&a_{3}\leftarrow\neg a_{1}\wedge\neg a_{2}\\ b_{1}\leftarrow\neg b_{2}\wedge\neg b_{3},&b_{2}\leftarrow\neg b_{3}\wedge\neg b_{1},&b_{3}\leftarrow\neg b_{1}\wedge\neg b_{2}\\ c_{1}\leftarrow\neg c_{2}\wedge\neg c_{3},&c_{2}\leftarrow\neg c_{3}\wedge\neg c_{1},&c_{3}\leftarrow\neg c_{1}\wedge\neg c_{2}\\ d_{1}\leftarrow\neg d_{2}\wedge\neg c_{3},&d_{2}\leftarrow\neg d_{3}\wedge\neg d_{1},&d_{3}\leftarrow\neg d_{1}\wedge\neg d_{2}\end{array}\right.

(45)

•

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

\displaystyle C_{1}

\displaystyle=

\displaystyle\left\{\;\begin{array}[]{lllll}\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}\right.

(51)

Assuming an ordering of atoms $\{a_{1},a_{2},a_{3},\ldots,d_{1},d_{2},d_{3}\}$ , the normal logic program $P_{1}$ shown in (45) is matricized to $P_{1}=(D_{1},Q_{1})$ where $D_{1}$ is a $(12\times 12)$ binary identity matrix (because there are 12 atoms and each atom has just one rule) and $Q_{1}$ is a $(12\times 24)$ binary matrix shown in (54). Constraints listed in (51) are a matricized to a ( $15\times 12$ ) constraint matrix $Q_{C_{1}}$ (57). In (54) and (57), $a$ for example stands for a triple $(a_{1}\,a_{2}\,a_{3})$ and $\neg a$ for $(\neg a_{1}\,\neg a_{2}\,\neg a_{3})$ .

\displaystyle Q_{1}\;

\displaystyle=

\displaystyle\begin{array}[]{ll}\begin{matrix}\;\;\;\;a&\;b&\;c&\;d&\hskip 1.00006pt\neg a&\hskip 1.00006pt\neg b&\hskip 1.00006pt\neg c&\hskip 1.00006pt\neg d\\ \vspace{-0.8em}\end{matrix}\\ \begin{bmatrix}&&&&\hskip 30.00005ptH_{3}&&&&\\ &&&&&H_{3}&&&\\ &&&&&&H_{3}&&\\ &&&&&&&H_{3}&\end{bmatrix}&\;\;\mbox{where}\;\;H_{3}\;=\;\begin{bmatrix}0&1&1\\ 1&0&1\\ 1&1&0\\ \end{bmatrix}\end{array}

(54)

\displaystyle Q_{C_{1}}\;

\displaystyle=

\displaystyle\begin{array}[]{ll}\begin{matrix}\;\;\;\;\;a&\hskip 6.00006ptb&\hskip 6.00006ptc&\hskip 5.0ptd\\ \vspace{-0.8em}\end{matrix}\\ \begin{bmatrix}E_{3}&E_{3}&&\\ E_{3}&&E_{3}&\\ &E_{3}&E_{3}&\\ &E_{3}&&E_{3}\\ &&E_{3}&E_{3}\\ \end{bmatrix}&\;\;\mbox{where}\;\;E_{3}\;=\;\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\\ \end{bmatrix}\end{array}

(57)

We run Algorithm 1 on program $P_{1}$ with constraints $C_{1}$ to find a supported model (solution) of $P_{1}$ satisfying $C_{1}$ ¹³¹³13 Since $P_{1}$ is a tight program, every supported model of $P_{1}$ is a stable model and vice versa. .

Table 1: Time and the number of solutions

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

To measure time to find a model, we conduct ten trials¹⁴¹⁴14 One trial consists of $max\_itr\times max\_try$ parameter updates. of running Algorithm 1 with $max\_try=20$ , $max\_itr=50$ , $\ell_{2}=\ell_{3}=0.1$ and take the average. The result is 0.104s(0.070) on average. Also to check the ability of finding different solutions, we perform ten trials of Algorithm 1¹⁵¹⁵15 without another solution constraint introduced in Section 5.2 and count the number of different solutions in the returned solutions. #solutions in Table 1 is the average of ten such measurements. Due to naive implementation, computation is slow but the number of different solutions, 5.2 on average, seems rather high considering there are six solutions.

Next we check the scalability of our approach by a simple problem. We consider the 3-coloring of a cycle graph like (a) in Fig. 2. In general, given a cycle graph that has $n$ nodes, we encode its 3-coloring problem as in the previous example by a matricized program $P_{2}=(D_{2},Q_{2})$ and a constraint matrix $Q_{C_{2}}$ where $D_{2}(3n\times 3n)$ is an identity matrix and $Q(3n\times 6n)$ and $Q_{C_{2}}(3n\times 6n)$ represent respectively rules and constraints. There are $2^{n}+2(-1)^{n}$ solutions ( $n\geq 3$ ) in $2^{3n}$ possible assignments for $3n$ atoms. So the problem will be exponentially difficult as $n$ goes up.

The graph (b) in Fig. 2 is an example of convergence curve of $J_{SU+c}$ by Algorithm 1 with $n=10,max\_try=100,max\_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\_itr=50$ iterations of update of continuous assignment vector ${\mathbf{u}}$ . Then perturbation is given to ${\mathbf{u}}$ which causes a small jump of $J_{SU+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{u}}$ to a binary vector ${\mathbf{u}}_{*}$ .

The graph (c) in Fig. 2 shows the scalability of computation time to find a solution up to $n=10000$ . We set $max\_try=100,max\_itr=2000$ and plot the average of ten measurements of time to find a solution. The graph seems to indicate good linearity w.r.t. $n$ up to $n=10000$ .

5.2 The Hamiltonian cycle problem, precomputation and another solution constraint

A Hamiltonian cycle (HC) is a cycle in a graph that visits every vertex exactly once and the Hamiltonian cycle 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 [19]. Later a way of encoding by a program that is not tight but tight on its completion models is proposed [13]. We here introduce yet another encoding by a tight ground program inspired by SAT encoding proposed in [30] where Zhou showed that the problem is solvable by translating six conditions listed in Fig. 3 into a SAT problem¹⁶¹⁶16 Actually he listed seven conditions to be encoded as a SAT problem. However, one of them is found to be redundant and we use the remaining six conditions. .

	conditions		meaning
(1)	${\rm one\_of}(H_{i,j_{1}},\ldots,H_{i,j_{k}})$	:	one of outgoing edges $\{i\rightarrow j_{1},\ldots,i\rightarrow j_{k}\}$ from vertex $i$ is in an HC
(2)	$U_{j,q}\leftarrow H_{i,j}\wedge U_{i,q-1}$	:	if edge $i\rightarrow j$ is in an HC and vertex $i$ is visited at time $q-1$ ,
			vertex $j$ is visited at time $q$ $(1\leq i,j,q\leq N)$
(3)	$U_{1,1}$	:	vertex 1 is visited at time 1 (starting point)
(4)	${\rm one\_of}(H_{i_{1},j},\ldots,H_{i_{k},j})$	:	one of incoming edges $\{i_{1}\rightarrow j,\ldots,i_{k}\rightarrow j\}$ to vertex $j$ is in an HC
(5)	$\leftarrow H_{i,1}\wedge\neg U_{i,N}$	:	if edge $i\rightarrow 1$ is in an HC, vertex $i$ is visited at time $N$ $(2\leq i\leq N)$
(6)	${\rm one\_of}(U_{i,1},\ldots,U_{i,N})$	:	vertex $i$ is visited once ( $1\leq i\leq N$ )

Figure 3: Conditions for SAT encoding

In what follows, we assume vertices are numbered from $1$ to $N$ = the number of vertices in a graph. We use $i\rightarrow 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\leq j,q\leq N$ ) and ${\rm one\_of}(a_{1},\ldots,a_{k})$ means that one of $\{a_{1},\ldots,a_{k}\}$ is true. We translate these conditions into a program $P_{3}=\{(1),(2),(3)\}$ and constraints $C_{3}=\{(4),(5),(6)\}$ . To be more precise, the first condition $(1)$ is translated into a tight program just like a program $P_{1}$ (45). The conditions $\{(2),(3)\}$ constitute a tight definite program. Constraints $C_{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},\ldots,L_{k}$ are literals. A set of $U_{j,q}$ atoms contained in a stable model of $P_{2}$ satisfying $C_{2}$ gives an HC.

We apply the above encoding to a simple Hamiltonian cycle problem for a graph $G_{2}$ in Fig. 4¹⁷¹⁷17 $G_{2}$ is taken from: Section 4.2 in A User’s Guide to gringo, clasp, clingo, and iclingo (http://wp.doc.ic.ac.uk/arusso/wp-content/uploads/sites/47/2015/01/clingo_guide.pd). . There are six vertices and six HCs¹⁸¹⁸18 They are $1\rightarrow 2\rightarrow 5\rightarrow 6\rightarrow 3\rightarrow 4\rightarrow 1$ , $1\rightarrow 2\rightarrow 6\rightarrow 3\rightarrow 5\rightarrow 4\rightarrow 1$ , $1\rightarrow 2\rightarrow 6\rightarrow 5\rightarrow 3\rightarrow 4\rightarrow 1$ , $1\rightarrow 3\rightarrow 5\rightarrow 6\rightarrow 2\rightarrow 4\rightarrow 1$ , $1\rightarrow 4\rightarrow 2\rightarrow 5\rightarrow 6\rightarrow 3\rightarrow 1$ , $1\rightarrow 4\rightarrow 2\rightarrow 6\rightarrow 5\rightarrow 3\rightarrow 1$ . . To solve this HC problem, we matricize $P_{3}$ and $C_{3}$ . There are 36 $H_{i,j}$ atoms ( $1\leq i,j\leq 6$ ) and 36 $U_{j,q}$ atoms ( $1\leq j,q\leq 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 $(D_{3},Q_{3})$ where $D_{3}$ is a $72\times 197$ binary matrix for disjunctions¹⁹¹⁹19 For example, for each $U_{j,q}$ ( $1\leq j,q\leq 6$ ), condition (2) generates six rules $\{U_{j,q}\leftarrow H_{i,j}\wedge U_{i,q-1}\mid 1\leq i\leq 6\}$ . and $Q_{3}$ is a $197\times 144$ matrix for conjunctions (rule bodies). Likewise $C_{3}=\{(4),(5),(6)\}$ is translated into a constraint matrix $Q_{C_{3}}$ which is a $67\times 144$ binary matrix. Then our task is to find a root ${\mathbf{u}}$ of $J_{SU+c}$ (3.5) constructed from these $D_{3}$ , $Q_{3}$ and $Q_{C_{3}}$ in a 72 dimensional vector space by minimizing $J_{SU+c}$ to zero.

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

We run Algorithm 1 on $P_{3}=(D_{3},Q_{3})$ with $Q_{C_{3}}$ (no precomputation) and also on $P^{\prime}_{3}=(D^{\prime}_{3},Q^{\prime}_{3})$ with $Q^{\prime}_{C_{3}}$ (precomputation) using $max\_try=20$ , $max\_itr=200$ and $\ell_{2}=\ell_{3}=0.1$ and measure time to find a solution, i.e. stable model satisfying constraints. The result is shown by Table (b) in Fig. 4 as time(s) where time(s) is an average of ten trials. Figures in the table, 2.08s vs. 0.66s, 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\ldots\wedge L_{72}$ of literals, we add a new constraint $\leftarrow L_{1}\wedge\ldots\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 ten 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, 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, $LF_{max}$ and $LF_{min}$ , introduced in the previous section work. We use an artificial non-tight program $P_{4\_n}$ (with no constraint) shown below that have exponentially many loops [12].

\displaystyle P_{4\_n}

\displaystyle=

\displaystyle\begin{array}[]{ll}\left\{\;\begin{array}[]{lllll}a_{0}&\leftarrow&a_{1}\wedge\ldots\wedge a_{n}\\ a_{0}&\leftarrow&\neg a_{n+1}\\ \ldots\\ a_{2i-1}&\leftarrow&a_{0}\vee a_{2i}&\mbox{for $i:1\leq i\leq n/2$}\\ a_{2i}&\leftarrow&a_{0}\vee a_{2i-1}&\mbox{for $i:1\leq i\leq n/2$}\\ \ldots\\ a_{n+1}&\leftarrow&a_{n+1}\\ \end{array}\right.\end{array}

For an even $n$ , $P_{4\_n}$ program has $n+2$ atoms $\{a_{0},a_{1},\ldots,a_{n},a_{n+1}\}$ , $2^{n/2}+1$ supported models and one stable model $\{a_{0},a_{1},\ldots,a_{n}\}$ . There are $n/2+1$ minimal loops $\{a_{1},a_{2}\},\ldots,$ $\{a_{n-1},a_{n}\},$ $\{a_{n+1}\}$ and a maximal loop $\{a_{0},a_{1},\ldots,a_{n}\}$ . The set of loop formulas for LF heuristics are computed as follows.

	$\displaystyle LF_{max}$	$\displaystyle=$	$\displaystyle\{(a_{0}\wedge a_{1}\wedge\ldots\wedge a_{n})\rightarrow\neg a_{n+1},\,a_{n+1}\rightarrow\perp\}$
	$\displaystyle LF_{max}$	$\displaystyle=$	$\displaystyle\{(a_{1}\wedge a_{2})\rightarrow a_{0},\ldots,(a_{n-1}\wedge a_{n})\rightarrow a_{0},\,a_{n+1}\rightarrow\perp\}$

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

After translating $P_{4\_n}$ into a matricized program $P_{4\_n}=(Q_{4\_n},D_{4\_n})$ where $Q_{4\_n}$ is a $(2n+3)\times(2n+4)$ binary matrix and $D_{4\_n}$ is a $(n+2)\times(2n+3)$ binary matrix respectively, we compute a stable model of $P_{4\_n}$ for various $n$ by Algorithm 1 that minimizes $J_{SU+c+LF}$ (37) with coefficient $\ell_{3}=0$ for the constraint term (because of no use of constraints) using Jacobian $J_{a_{SU+c+LF}}$ (32).

Since all supported models of $P_{4\_n}$ except for one stable model are non-stable, even if $LF_{max}$ and $LF_{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 situation by the use of another solution constraint.

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)$

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 $LF_{max}$ or $LF_{min}$ heuristics, we repeatedly run Algorithm 1 with/without another solution constraint for $10^{4}$ trials with $n=4$ , $max\_try=20$ , $max\_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 ten times and take the average. The result is shown in Table 2.

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.09s 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,\ldots,50$ , we repeatedly run Algorithm 1 using $J_{SU+c+LF}$ with $max\_try=10,max\_itr=100$ on matricized $P_{4\_n}=(Q_{4\_n},D_{4\_n})$ (and no constraint matrix) to compute supported (stable) models. Coefficients in $J_{SU+c+LF}$ are set to $\ell_{2}=0.1,\ell_{3}=0,\ell_{4}=1$ . To be more precise, for each $n$ and each case of $LF_{max}$ , $LF_{min}$ , precomputation (without $\{LF_{max},LF_{min}\}$ ) and no $\{LF_{max},LF_{min},\mbox{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 ten times and take the average and obtain graphs in Fig. 5.

In Fig. 5, no_LF means no use of $\{LF_{max},LF_{min}\}$ heuristics. Also no_LF_pre means no_LF is applied to precomputed $P_{4\_n}$ ²⁰²⁰20 Precomputation takes 0.006s and removes only one stable false atom. So precomputation is not helpful in the current case. .

We can see from graph (a) in Fig. 5 that computation time is $LF_{min}>LF_{max}>\mbox{no\_LF}>\mbox{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 Fig. 5 tells us that $LF_{min}$ allows computing redundant non-stable models but the rest, $LF_{max}$ , no_LF and no_LF_pre, return a stable model without computing redundant non-stable models. This shows first that $LF_{max}$ works correctly to suppress the computation of non-stable models and second that the $LF_{min}$ heuristics works adversely, i.e. 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\_nk}$ . It is a modification of $P_{4\_n}$ by adding self-loops of $k$ atoms as illustrated by (a) in Fig. 6. The addition of self-loop causes the choice of $a_{n+j}$ ( $1\leq j\leq k$ ) being true or being false in the search process. $P_{5\_nk}$ has $(2^{n/2}-1)(2^{k}-1)+1$ supported models but has just one stable model $\{a_{0},a_{1},\ldots,a_{n}\}$ .

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

We run Algorithm 1 on $P^{\prime}_{5\_nk}$ with $\ell_{2}=\ell_{3}=0.1$ and $max\_try=10$ , $max\_itr=100$ at most 100 trials to measure time to find a stable model ten times for each $n=1000,\ldots,5000$ and take the average. At the same time, we also run Clingo (version 5.6.2) on $P_{5\_nk}$ and similarly measure time. Graph (b) in Fig. 6 is the result. It shows that as far as computing a stable model of $P_{5\_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},\ldots,a_{n+k}\}$ in the stable model of $P_{5\_nk}$ and Algorithm 1 run on the precomputed $P^{\prime}_{5\_nk}=(D^{\prime}_{5\_nk},Q^{\prime}_{5\_nk})$ detects the stable model only by thresholding ${\mathbf{u}}$ before starting any update of ${\mathbf{u}}$ . 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 [2] and [27]. 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\neg b,b\leftarrow\neg a\}$ . According to [2], the success rate w.r.t. $N$ of returning a supported model goes from one initially to almost zero at $N=64$ in [2] while it keeps one till $N=100$ in [27]. We tested the same program with $max\_try=20$ , $max\_itr=100$ and observed that the success rate keeps one till $N=10000$ .

Although our approach is non-probabilistic, i.e. purely linear algebraic, there are probabilistic differentiable approaches for ASP. Differentiable ASP/SAT [18] iteratively samples a stable model by an ASP solver a la ASSAT [14]. 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 [29] 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 [26] 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.

Sato and Kojima developed another differentiable approach which, without using ASP solvers, samples supported models of probabilistic normal logic programs [24]. They encode programs by matrices and realize the sampling of supported models in vector spaces by repeatedly computing a fixedpoint of some differentiable equations.

Independently of differentiable approaches mentioned so far, Tuan et al. adopted matrix encoding for propositional normal logic programs based on [22] and proposed to compute stable models in vector spaces by a generate-and-test manner using sparse representation [21].

Also though not computing stable models, Sato et al. developed a 2-stage approach for supported models. It first computes the least 3-valued model of a matricized program in vector spaces that contains undefined atoms, and then assigns true or false to them to obtain a supported model of the original program [25].

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 [14] 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.

Acknowledgments

This work is supported by JSPS KAKENHI Grant Number JP21H04905 and JST CREST Grant Number JPMJCR22D3..

References

[1] Apt, K.R., Blair, H.A., Walker, A.: Foundations of deductive databases and logic programming. chap. Towards a Theory of Declarative Knowledge, pp. 89–148 (1988)
[2] Aspis, Y., Broda, K., Russo, A., Lobo, J.: Stable and supported semantics in continuous vector spaces. In: Calvanese, D., Erdem, E., Thielscher, M. (eds.) Proceedings of the 17th International Conference on Principles of Knowledge Representation and Reasoning, KR 2020, Rhodes, Greece, September 12-18, 2020. pp. 59–68 (2020)
[3] Dowling, W.F., Gallier, J.H.: Linear-time Algorithms for Testing the Satisfiability of Propositional Horn Formula. Journal of Logic Programming 3, 267–284 (1984)
[4] Erdem, E., Lifschitz, V.: Tight Logic Programs. Theory and Practice of Logic Programming (TPLP) 3(4–5), 499–518 (2003)
[5] Fages, F.: Consistency of Clark’s completion and existence of stable models. Journal of Methods of Logic in Computer Science 1, 51––60 (1994)
[6] Ferraris, P., Lee, J., Lifschitz, V.: A generalization of the lin-zhao theorem. Annals of Mathematics and Artificial Intelligence 47, 79–101 (2006)
[7] Gebser, M., Kaufmann, B., AndréNeumann, Schaub, T.: clasp: A conflict-driven answer set solver. In: LPNMR. Lecture Notes in Computer Science, vol. 4483, pp. 260–265. Springer (2007)
[8] Gebser, M., Kaufmann, B., Schaub, T.: Conflict-driven answer set solving: From theory to practice. Artif. Intell. 187, 52–89 (2012)
[9] Gelfond, M., Lifshcitz, V.: The stable model semantics for logic programming pp. 1070–1080 (1988)
[10] Lierler, Y.: Cmodels: Sat-based disjunctive answer set solver. In: Proceedings of the 8th International Conference on Logic Programming and Nonmonotonic Reasoning. p. 447–451. LPNMR’05, Springer-Verlag, Berlin, Heidelberg (2005)
[11] Lifschitz, V.: What is answer set programming? In: Proceedings of the 23rd National Conference on Artificial Intelligence - Volume 3. p. 1594–1597. AAAI’08, AAAI Press (2008)
[12] Lifschitz, V., Razborov, A.: Why are there so many loop formulas? ACM Trans. Comput. Logic 7(2), 261–268 (apr 2006). https://doi.org/10.1145/1131313.1131316
[13] Lin, F., Zhao, J.: On tight logic programs and yet another translation from normal logic programs to propositional logic. In: Proceedings of the 18th International Joint Conference on Artificial Intelligence (IJCAI’03). p. 853–858. Morgan Kaufmann Publishers Inc. (2003)
[14] Lin, F., Zhao, Y.: Assat: computing answer sets of a logic program by sat solvers. Artificial Intelligence 157(1), 115–137 (2004)
[15] Liu, H., Wang, J.: A new way to enumerate cycles in graph. Advanced Int’l Conference on Telecommunications and Int’l Conference on Internet and Web Applications and Services (AICT-ICIW’06) pp. 57–59 (2006)
[16] Marek, V.W., Truszczyński, M.: Stable Models and an Alternative Logic Programming Paradigm, pp. 375–398. Springer Berlin Heidelberg (1999). https://doi.org/10.1007/978-3-642-60085-2_17
[17] Marek, W., V.S.Subrahmanian: The relationship between stable, supported, default and autoepistemic semantics for general logic programs. Theoretical Computer Science 103(2), 365–386 (1992)
[18] Nickles, M.: Differentiable SAT/ASP. In: Proceedings of the 5th International Workshop on Probabilistic Logic Programming, PLP 2018. pp. 62–74 (2018)
[19] Niemelä, I.: Logic programs with stable model semantics as a constraint programming paradigm 25(3–4), 241–273 (1999)
[20] Niemelä, I., Simons, P.: Smodels — an implementation of the stable model and well-founded semantics for normal logic programs. In: Dix, J., Furbach, U., Nerode, A. (eds.) Logic Programming And Nonmonotonic Reasoning. pp. 420–429. Springer Berlin Heidelberg (1997)
[21] Quoc, T.N., Inoue, K., Sakama, C.: Enhancing linear algebraic computation of logic programs using sparse representation. New Gen. Comput. 40(1), 225–254 (2022). https://doi.org/10.1007/s00354-021-00142-2
[22] Sakama, C., Inoue, K., Sato, T.: Linear Algebraic Characterization of Logic Programs. In: Proceedings of the 10th International Conference on Knowledge Science, Engineering and Management (KSEM2017), LNAI 10412, Springer-Verlag. pp. 520–533 (2017)
[23] Sarker, M.K., Zhou, L., Eberhart, A., Hitzler, P.: Neuro-symbolic artificial intelligence: Current trends (2021)
[24] Sato, T., Kojima, R.: In: Artificial Intelligence IJCAI 2019 International Workshops, Revised Selected Best Papers (Amal El Fallah Seghrouchni,David Sarne (Eds.)), Lecture Notes in Computer Science, vol. 12158, pp. 239–252. Springer (2020). https://doi.org/10.1007/978-3-030-56150-5
[25] Sato, T., Sakama, C., Inoue, K.: From 3-valued semantics to supported model computation for logic programs in vector spaces. In: Rocha, A.P., Steels, L., van den Herik, H.J. (eds.) Proceedings of the 12th International Conference on Agents and Artificial Intelligence, ICAART 2020, Volume 2, Valletta, Malta, February 22-24, 2020. pp. 758–765. SCITEPRESS (2020). https://doi.org/10.5220/0009093407580765
[26] Skryagin, A., Stammer, W., Ochs, D., Dhami, D.S., Kersting, K.: Neural-Probabilistic Answer Set Programming. In: Proceedings of the 19th International Conference on Principles of Knowledge Representation and Reasoning. pp. 463–473 (2022)
[27] Takemura, A., Inoue, K.: Gradient-based supported model computation in vector spaces. In: Gottlob, G., Inclezan, D., Maratea, M. (eds.) Logic Programming and Nonmonotonic Reasoning - 16th International Conference, LPNMR 2022, Genova, Italy, September 5-9, 2022, Proceedings. Lecture Notes in Computer Science, vol. 13416, pp. 336–349. Springer (2022)
[28] Tarjan, R.E.: Depth-first search and linear graph algorithms. SIAM Journal on Computing 1(2), 146–160 (1972)
[29] Yang, Z., Ishay, A., Lee, J.: Neurasp: Embracing neural networks into answer set programming. In: Bessiere, C. (ed.) Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI-20. pp. 1755–1762. International Joint Conferences on Artificial Intelligence Organization (2020)
[30] Zhou, N.F.: In pursuit of an efficient sat encoding for the hamiltonian cycle problem. In: proceeding of the Principles and Practice of Constraint Programming: 26th International Conference (CP 2020). pp. 585–602. Springer-Verlag (2020)