Enhancing Linear Algebraic Computation of
Logic Programs Using Sparse Representation

We consider a language $\mathscr{L}$ that contains a finite set of propositional variables. Given a logic program $P$, the set of all propositional variables appearing in $P$ is called the Herbrand base of $P$ (written $B_{\_}P$). A definite program is a finite set of rules of the form:

where $h$ and $b_{\_}i$ are propositional variables (atoms) in $\mathscr{L}$.

A rule $r$ is called an OR-rule if $r$ is the form:

where $h$ and $b_{\_}i$ are propositional variables in $\mathscr{L}$.

A standardized program is a finite set of rules that are either (1) or (2). Note that the rule (2) is a shorthand of $m$ rules: $h\leftarrow b_{\_}1$, $\ldots$, $h\leftarrow b_{\_}m$, so a standardized program is considered a definite program.

For each rule $r$ of the form (1) or (2), define $head(r)=h$ and $body(r)=\{b_{\_}1,\ldots,b_{\_}m\}$.¹¹1We assume $b_{\_}i\neq b_{\_}j\;\forall\;i\neq j$. A rule $r$ is called a fact if $body(r)=\emptyset$.

Any definite program can be transformed to a standardized program by introducing new propositional variables. That is, if there are several rules with the same head but different bodies: $h\leftarrow body(r_{\_}1)$, $\ldots$, $h\leftarrow body(r_{\_}k)$, then we replace all these rules by $h\leftarrow b_{\_}1\vee\ldots\vee b_{\_}k$, $b_{\_}1\leftarrow body(r_{\_}1)$,… $b_{\_}k\leftarrow body(r_{\_}k)$. In this paper, a program means a standardized program unless stated otherwise.

A set $I\subseteq B_{\_}P$ is an interpretation of $P$. An interpretation $I$ is a model of a standardized program $P$ if $\{b_{\_}1,\ldots,b_{\_}m\}\subseteq I$ implies $h\in I$ for every rule (1) in $P$, and $\{b_{\_}1,\ldots,b_{\_}m\}\cap I\neq\emptyset$ implies $h\in I$ for every rule (2) in $P$. A model $I$ is the least model of $P$ if $I\subseteq J$ for any model $J$ of $P$. A mapping $T_{\_}P:\,2^{B_{\_}P}\rightarrow 2^{B_{\_}P}$ (called a $T_{\_}P$-operator) is defined as: $T_{\_}P(I)=\{\,h\,\mid\,h\leftarrow b_{\_}1\wedge\cdots\wedge b_{\_}m\in P\;\mbox{and}\;\{b_{\_}1,\ldots,b_{\_}m\}\subseteq I\,\}\;\cup\;\{\,h\,\mid\,h\leftarrow b_{\_}1\vee\cdots\vee b_{\_}n\in P\;\mbox{and}\;\{b_{\_}1,\ldots,b_{\_}n\}\cap I\neq\emptyset\,\}.$

The powers of $T_{\_}P$ are defined as: $T_{\_}P^{k+1}(I)=T_{\_}P(T_{\_}P^{k}(I))$ $(k\geq 0)$ and $T_{\_}P^{0}(I)=I$. Given $I\subseteq B_{\_}P$, there is a fixpoint $T_{\_}P^{n+1}(I)=T_{\_}P^{n}(I)$ $(n\geq 0)$. For a definite program $P$, the fixpoint $T_{\_}P^{n}(\emptyset)$ coincides with the least model of $P$ [18].

1. 1.

   $a_{\_}{ij_{\_}k}=\frac{1}{m}\;\;(1\leq k\leq m;\,1\leq i,j_{\_}k\leq n)$ if  $p_{\_}i\leftarrow p_{\_}{j_{\_}1}\wedge\cdots\wedge p_{\_}{j_{\_}m}$ is in $P$;

2. 2.

   $a_{\_}{ij_{\_}k}=1\;\;(1\leq k\leq l;\,1\leq i,j_{\_}k\leq n)$ if  $p_{\_}i\leftarrow p_{\_}{j_{\_}1}\vee\cdots\vee p_{\_}{j_{\_}l}$ is in $P$;

3. 3.

   $a_{\_}{ii}=1$ if $p_{\_}i\leftarrow$ is in $P$;

4. 4.

   $a_{\_}{ij}=0$, otherwise.

$M_{\_}P$ is called a program matrix. We write ${\sf row}_{\_}i(M_{\_}P)=p_{\_}i$ and ${\sf col}_{\_}j(M_{\_}P)=p_{\_}j$ $(1\leq i,j\leq n)$.

To better understand Definition 2, let’s consider a concrete example.

Sakama et al. further define representation of interpretation by using interpretation vector (Definition 3). This vector is used to store the truth value of all propositions in $P$. The starting point of interpretation vector is defined as the initial vector (Definition 4).

In order to compute the least model in vector space, Sakama et al. proposed an algorithm which is equivalent to the result of computing least models by the $T_{\_}P$-operator. This algorithm is presented in Algorithm 1.

Normal programs can be transformed to definite programs as introduced in [2]. Therefore, we transform normal programs to definite programs before encoding them in matrices.

$P$ is transformed to a definite program by rewriting the above rule to the following form:

where $\overline{b}_{\_}i$ is a new proposition associated with $b_{\_}i$.

In this part, we denote $P$ as a normal program with an interpretation $I\subseteq B_{\_}P$. The positive form $P^{+}$ of $P$ is obtained by applying the above transformation. Since a definite program $P^{+}$ is transformed to its standardized program, then we can apply Algorithm 1 to compute the least model. [2] proved that if $P$ is a normal program, $I$ is a stable model of $P$ iff $I^{+}$ is the least model of $P^{+}\cup\overline{I}$. We should note that $I^{+}$ is interpretation of $P^{+}$ which is a definite program. We can obtain $P^{+}$ by applying Algorithm 1 to the transformed program $P^{+}$. Define $\bar{I}=\{\bar{p}\ |\ p\in B_{\_}P\setminus I\}$, then $I^{+}=I\cup\bar{I}$.

1. 1.

   $a_{\_}{ii}=1$ for $n+1\leq i\leq m$;

2. 2.

   $a_{\_}{ij}=0$ for $n+1\leq i\leq m$ and $1\leq j\leq m$ such that $i\neq j$;

3. 3.

   Otherwise, $a_{\_}{ij}$ ($1\leq i\leq n$; $1\leq j\leq m$) is encoded as in Definition 2.

$M_{\_}P$ is called a program matrix. We write ${\sf row}_{\_}i(M_{\_}P)=p_{\_}i$ and ${\sf col}_{\_}j(M_{\_}P)=p_{\_}j$ $(1\leq i,j\leq n)$.

Instead of the initial vector in the case of definite programs, the notion of an initial matrix is introduced to encode positive and negative facts in a program.

1. 1.

   each row of $M_{\_}0$ corresponds to each element of $B_{\_}P$ in a way that $row_{\_}i(M_{\_}0)=p_{\_}i$ for $1\leq i\leq n$ and $row_{\_}i(M_{\_}0)=\overline{q}_{\_}i$ for $n+1\leq i\leq m$;

2. 2.

   $a_{\_}{ij}=1$ ($1\leq i\leq n$, $1\leq j\leq h$) iff a fact $q_{\_}i\leftarrow$ is in $P$; otherwise $a_{\_}{ij}=0$;

3. 3.

   $a_{\_}{ij}=0$ ($n+1\leq i\leq m$, $1\leq j\leq h$) iff a fact $q_{\_}i$ is in $P$; otherwise, there are two possibilities $0$ and $1$ for $a_{\_}{ij}$, so it is either $0$ or $1$.

Each column of $M_{\_}0$ is a potential stable model in the first stage. We update $M_{\_}0$ by applying matrix multiplication with the matrix representation obtained by Definition 7 as $M_{\_}{k+1}=\theta(M_{\_}PM_{\_}k)$. Then, the algorithm for computing the stable models is presented in Algorithm 2.

This method requires extra steps on transforming and finding stable models of a program. In addition, the initial matrix size grows exponentially by the number of negations $m-n$. Therefore this representation requires a lot of memory and the algorithm performs considerably slower than the method for definite programs if there are many negations appear in the program. But we will later show that this method still has the advantage when there are a small number of negations.

A sparse matrix is a matrix in which most of the elements are zero. The level of sparseness is measured by sparsity which equals the number of zero-valued elements divided by the total number of elements [5]. Because there are a large number of zero elements in sparse matrices, we can save the computation by ignoring these zero values [9]. According to the conversion method of linear algebraic approach, we can calculate the sparsity of a program $P$ ²²2We only consider the programs in Definition 2 and Definition 7.. This calculation is done by counting the number of nonzero-valued elements of each rule in $P$, then let $1$ minus the fraction of the number of nonzero-valued elements and the matrix size.

By definition, the sparsity of a program $P$ is computed by the following equation:

where $n$ is the number of elements in $B_{\_}P$ and $|B_{\_}r|$ is the length of body of rule $r$.

Accordingly, the representation matrix becomes a high level of sparsity if matrix size becomes larger while the length of body rule is insignificant. In fact, a rule $r$ in a logic program rarely has a body length approx $n$, therefore, $|B_{\_}r|\ll n$. In short, we can say that the matrix representation of a logic program according to the linear algebraic approach is highly sparse.

There are several sparse matrix representations. Compressed Sparse Row (CSR) is one of the most efficient, robust and widely adopted by many programming libraries among those [5]. Hence, in this paper, we propose CSR for representing logic programs.

In order to understand the CSR format, firstly we need to mention the Coordinate (COO) format which is simple using 2 arrays of coordinates and 1 array of values. The length of these arrays is equal to the number of nonzero elements. The first array stores the row index of each value, and the second array stores the row and column indices of each value, while the third array stores the values in the original matrix. We can imagine that the $i^{th}$ nonzero element in a matrix is represented by a 3-tuple extracted from these 3 arrays at index $i$. Example 3 illustrates sparse representation in COO format for the program $P$ in Example 1. We should note that in Example 3, zero-based indexing is used.

Noticeably, in the row index array, it is possible for a value to be repeated consecutively because the nonzero elements may appear in the same row for many times. We may reduce the size of the row index array by considering CSR format. In this format, while the column index and the value arrays remain the same, we compress the row index array by storing the index of the row only where nonzero elements appear. That means we do not need to store 2 consecutive $0$ and 2 consecutive $5$ as in Example 3. Instead, we store the index of the next row, then finally point the last index to the end of the row (which equals to the number of nonzero elements). Concretely in the row index array, the first element is starting index which is 0. The last element is an extra element to indicate the end of this array which is equal to the number of nonzero elements. We need 2 consecutive values in row index array to extract the nonzero elements in this row. To be specific, we need to interpret $row\_start$ and $row\_end$ of the $i^{th}$ row from the compressed value in $row\_index$ array: $row\_start_{\_}i=row\_index[i],row\_end_{\_}i=row\_index[i+1]$.

Example 4 illustrates this method. For the first row ($i=0$), we have $row\_start_{\_}0=0,row\_end_{\_}i=2$, then we extract 2 values $0$ and $1$ for the nonzero element in the first row. These $start$ and $end$ will be used to extract column index and value of nonzero elements. Similarly, the second row ($i=1$), we have $row\_start_{\_}1=2,row\_end_{\_}1=3$ then we have only one nonzero element at index 2. Continue this interpretation until we reach the final row ($i=6$), we have $row\_start_{\_}6=8,row\_end_{\_}6=10$ then we extract last two nonzero elements at index 8 and 9 for the final row.

As we can see in Example 4, the row index array now has only 8 indices rather than 10 in Example 3. We save storing repeatedly indices in the row index array by storing only the position where it starts and ends. This phenomenon does not always encourage the reduction in terms of array size. Imagine the situation where each row has only one nonzero element, then we save nothing with this representation. Actually for a sparse matrix of the size $m\times n$, the CSR format saves on memory only when $nnz<(m(n-1)-1)/2$ (where $nnz$ is number of nonzero elements). Fortunately, in case of linear algebraic method, each rule normally has many atoms in the body, therefore the matrix representation has many nonzero elements in a single row. That is why CSR format could considerably save memory. In fact, we can save up to $20\%$ of the size of the row index array using CSR format. Accordingly, in our method, we propose to implement CSR rather than COO not only because it saves more memory, but also because it is widely adopted by many programming libraries.

Because a logic program $P$ is highly sparse, applying Algorithm 1 and Algorithm 2 on sparse representation is remarkably faster than the dense matrix. Moreover, sparse representation saves the memory space as well, therefore enabling the ability to deal with large scale KBs. We are going to reveal the performance gain by using sparse representation in the next section.

Note that only the matrix representation of the program is sparse, while the initial vector and the initial matrix are not sparse. Thus, in our methods, we will keep the dense format for interpretation matrices.

The final results on definite programs are illustrated in Table 2. We can see in the results, Dense matrix method is the slowest method and being unable to run with very large programs. Overall, Sparse matrix method is very efficient which is $10\sim 15$ faster than Clasp. We should mention that all the codes are executed on single-threaded CPU without using GPU boost or any other parallel computing techniques.

The benchmark on denser matrix are presented in Table 3. As can be seen in the results, denser matrices require more computation for Sparse matrix method while it does not affect the same scale on other competitors. Despite of that fact, Sparse matrix method still holds the first place in this benchmark. In terms of analyzing the sparseness level of logic programs, we hardly find a program in which the sparsity is less than 0.97. This observation strongly encourages the use of sparse representation for logic programs.

We next show the comparison for computing transitive closure. We assume that a dataset contains $edges$ (tuples of $nodes$), then first perform grounding two rules of defining $path$. The obtained results are demonstrated in Table 4. In this non-randomized problem, we can see that the matrix representations are very sprase. Therefore, it is no doubt that Sparse matrix method outperforms Dense matrix method. Accordingly, we only highlight the efficiency of sparse representation and omit the dense matrix approach. Surpringly, Sparse matrix method surpasses Clasp once again in this experiment by a large margin.

As can be witnessed in the results, Dense matrix method is the slowest, even slower than $T_{\_}p$-operator, in terms of computation time due to wasting of computation on a huge amount of zero elements. This could be explained by the high level of sparsity of logic programs provided in Table 2, Table 3 and Table 4. Moreover, large dense matrices consume a huge amount of memory therefore the method is unable to run with large scale matrix size. Overall, sparse matrix method is effective in computing the fixpoint of definite programs. On the other hand, the performance would be improved if we use GPU accelerated code and exploit parallel computing power. The results indicates a potential for logical inference using an algebraic method.

In our current method, the number of columns in the initial matrix (Definition 8) grows exponentially by the number of negations, we limit the number of negations in this benchmark by $8$ as specified in the experiment setup.

First, we perform benchmarks on normal programs which has $0.99$ sparsity level. Table 5 illustrates the execution time in detail. As can be witnessed in the results, Sparse matrix method is still faster than Clasp but with a smaller scale it did in definite programs. It is needed to mention that the initial matrix size is remarkably larger due to the limitation of representation. We have to initialize all possible combinations of an atom which appeared with its negation form in the program. There is no doubt that with a larger number of negations, the space complexity of linear algebraic method is exponential. Accordingly, the performance of Sparse matrix method is only better than Clasp with a small fraction of negations.

In the next experiments, we compare different methods on denser matrix. Table 6 presents the data for this benchmark. Once again, with a limited number of negations, Sparse matrix method holds the winner position.

Noticeably, execution time on normal programs is generally greater than that on definite programs. This is obvious because we have a larger size of initial matrices as well as the need of extra computation on transforming and finding the least models as described in Algorithm 2. Then the weakness of the linear algebraic method is that we have to deal with all combinations of truth assignments in order to compute the stable model. Accordingly, the column size of the initial matrix exponentially increases by the number of negations. Thus, in the benchmark on randomized programs, we limit the number of negations for all benchmarks so that the matrix can fit in memory. This limitation will become clearer in real problems which have many negations. This is a major problem that we are investigating to do further research.